COVID-19 in Louisiana - A Time Series Analysis
Dustin Bracy
12/5/2020
Introduction
My name is Dustin Bracy, and this is a time series analysis of COVID-19 data in Louisiana and the United States spanning an 8 month period from the declaration of the pandemic on March 11, 2020 through November 11, 2020. Below I will explore the difference between new daily positive cases and positivity rates, and why one may be more advantageous. Next I’ll walk through univiariate and multivariate analyses of both Louisiana and national numbers. Finally I’ll wrap up with a comparison of model performance for each use case.
Acquire and Load the Data
The COVID-19 data used here comes from two primary sources: the Louisiana Department of Health Coronavirus website, and The Covid Tracking Project website. Additional temperature, windspeed and other related data for Louisiana was collected manually from historical data for the MSY (New Orleans regional airport) section of Weather Underground.
The data comes in pretty good shape, with only a handful of missing values in the US dataset. Fortunately, there are no missing values in response variables of interest (i.e. the total number of daily tests and the number of daily positive cases).
The data ingestion pipeline is fairly simple: download the data directly from the api or website, remove data outside of the March 11 to November 11 window, and add positivity metrics.
# Get latest US Data & sort by date ascending
download.file('https://covidtracking.com/data/download/national-history.csv',
destfile = './data/national-history.csv')
us <- read_csv('./data/national-history.csv')##
## ── Column specification ─────────────────────────────────────────────────────
## cols(
## date = col_date(format = ""),
## death = col_double(),
## deathIncrease = col_double(),
## inIcuCumulative = col_double(),
## inIcuCurrently = col_double(),
## hospitalizedIncrease = col_double(),
## hospitalizedCurrently = col_double(),
## hospitalizedCumulative = col_double(),
## negative = col_double(),
## negativeIncrease = col_double(),
## onVentilatorCumulative = col_double(),
## onVentilatorCurrently = col_double(),
## positive = col_double(),
## positiveIncrease = col_double(),
## recovered = col_double(),
## states = col_double(),
## totalTestResults = col_double(),
## totalTestResultsIncrease = col_double()
## )us <- us[order(us$date),]
# Get latest LA Data & consolidate to state level
download.file('https://ldh.la.gov/assets/oph/Coronavirus/data/LA_COVID_TESTBYDAY_PARISH_PUBLICUSE.xlsx',
destfile = './data/LA_COVID_TESTBYDAY_PARISH_PUBLICUSE.xlsx')
la <- readxl::read_xlsx('./data/LA_COVID_TESTBYDAY_PARISH_PUBLICUSE.xlsx')
la <- la %>%
group_by(date = `Lab Collection Date`) %>%
summarise(
tests = sum(`Daily Test Count`),
positive = sum(`Daily Positive Test Count`),
negative = sum(`Daily Negative Test Count`),
cases = sum(`Daily Case Count`)
)## `summarise()` ungrouping output (override with `.groups` argument)# Create positivity features
la$positivity = la$positive / la$tests
us$positivity <- us$positiveIncrease / us$totalTestResultsIncrease
# Get temp data from WeatherUnderground.com (historical data for MSY)
## Note Nov8 data missing from site, so the imputed average of Nov7+9 was taken
xdf <- read_csv('./data/temp.csv')##
## ── Column specification ─────────────────────────────────────────────────────
## cols(
## Date = col_character(),
## MaxTempF = col_double(),
## AvgTempF = col_double(),
## MinTempF = col_double(),
## MaxDewPointF = col_double(),
## AvgDewPointF = col_double(),
## MinDewPointF = col_double(),
## MaxHumidityPercent = col_double(),
## AvgHumidityPercent = col_double(),
## MinHumidityPercent = col_double(),
## MaxWindSpeedMPH = col_double(),
## AvgWindSpeedMPH = col_double(),
## MinWindSpeedMPH = col_double(),
## MaxPressureHg = col_double(),
## AvgPressureHg = col_double(),
## MinPressureHg = col_double(),
## PrecipitationIN = col_double(),
## Weekday = col_character(),
## Gathering = col_character()
## )xdf$GatheringBinary <- as.numeric(ifelse(is.na(xdf$Gathering), 0,1))
# Trim all datasets to span March 11, 2020 through Nov 11, 2020
la <- la[11:256,]
us <- us[50:295,]
xdf <- xdf[11:256,]
# Check for missing values
la %>% is.na() %>% summary()## date tests positive negative
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:246 FALSE:246 FALSE:246 FALSE:246
## cases positivity
## Mode :logical Mode :logical
## FALSE:246 FALSE:246us %>% is.na() %>% summary()## date death deathIncrease inIcuCumulative
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:246 FALSE:246 FALSE:246 FALSE:232
## TRUE :14
## inIcuCurrently hospitalizedIncrease hospitalizedCurrently
## Mode :logical Mode :logical Mode :logical
## FALSE:231 FALSE:246 FALSE:240
## TRUE :15 TRUE :6
## hospitalizedCumulative negative negativeIncrease onVentilatorCumulative
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:246 FALSE:246 FALSE:246 FALSE:225
## TRUE :21
## onVentilatorCurrently positive positiveIncrease recovered
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:232 FALSE:246 FALSE:246 FALSE:232
## TRUE :14 TRUE :14
## states totalTestResults totalTestResultsIncrease positivity
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:246 FALSE:246 FALSE:246 FALSE:246
## xdf %>% dplyr::select(!Gathering) %>% is.na() %>% summary()## Date MaxTempF AvgTempF MinTempF
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:246 FALSE:246 FALSE:246 FALSE:246
## MaxDewPointF AvgDewPointF MinDewPointF MaxHumidityPercent
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:246 FALSE:246 FALSE:246 FALSE:246
## AvgHumidityPercent MinHumidityPercent MaxWindSpeedMPH AvgWindSpeedMPH
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:246 FALSE:246 FALSE:246 FALSE:246
## MinWindSpeedMPH MaxPressureHg AvgPressureHg MinPressureHg
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:246 FALSE:246 FALSE:246 FALSE:246
## PrecipitationIN Weekday GatheringBinary
## Mode :logical Mode :logical Mode :logical
## FALSE:246 FALSE:246 FALSE:246Comparing Positivity Rate vs New Positive Cases
Positivity rate - According to John Hopkins University, the positivity rate is a new measure developed to identify either trending increases in transmission rates, or overall drops in testing rates, which can indicate a shortage of tests. It measures the ratio of diagnostic tests which are actually positive (i.e. # of positive tests / # of total tests).
The idea behind this test is that a high positivity score should give leaders indication that either infection is increasing at an unusual rate, or not enough tests are available. The World Health Organization has recommended nations adopt a rule of thumb to keep populations under approximately 5% positivity rate before considering re-opening or relaxing social restrictions.
The positivity rate improves upon strictly measuring the total number of cases or positive results by better accounting for new tests. As time goes on, more individuals will be re-tested, including high risk individuals such as health care workers, who may receive tests and positive results more frequently than others. The positivity rate accounts for this by measuring the total number of positive cases against tests, which will show less bias than strictly reporting on positive cases, which can be inflated by duplicate tests and non-unique cases. A better way to control this bias is to measure the people testing positive over the people who have been tested, although this data is much more difficult to collect, due to HIPAA and privacy concerns with linking data to individuals.
Sources https://www.jhsph.edu/covid-19/articles/covid-19-testing-understanding-the-percent-positive.html https://www.mprnews.org/story/2020/08/12/behind-the-numbers-what-does-a-covid19-positivity-rate-really-mean
Below is a plot that shows the differences between new positive cases vs the positivity rate. You can start to see the explosive nature of the positive case metric. I determined the ratio of Louisiana’s population compared to the US population by using census.gov population values from Jul 1, 2019.
#census.gov reported populations on Jul 1, 2019: Louisiana / USA
lapopratio <- 4648794 / 328239523
positivity_df <- data.frame(cbind(la[c('date', 'positivity')], us['positivity']))
names(positivity_df) = c('date','LA','US')
positive_df <- data.frame(cbind(la[c('date', 'positive')], us['positiveIncrease']))
names(positive_df) = c('date','LA','US')
positive_df$US = positive_df$US * lapopratio
grid.arrange(
gather(positivity_df, key = Region, value = 'Positivity Rate', -date) %>%
ggplot(aes(date, `Positivity Rate`)) +
geom_line(aes(color=Region), size=.75) +
labs(title='Louisiana vs US Positivity Rate',
x = 'Date',
y='Positivity Rate'),
gather(positive_df, key = Region, value = 'Positive Tests', -date) %>%
ggplot(aes(date, `Positive Tests`)) +
geom_line(aes(color=Region), size=.75) +
labs(title='Louisiana vs US Daily New Positive Tests',
x = 'Date',
y='Positive Tests',
subtitle = '*US numbers scaled to 1.42% based on LA population density'),
ncol = 2)
Evaluation Methods
To compare the prediction performance of our models to reality, we want to train a model and hold back a test set of a given period of days that I’ll call a horizon. We will sum the squared difference between each predicted value and actual value, and then take the average to generate the Average Squared Error (ASE). This metric will be the base comparison for all models, and generally the model with the lowest ASE score is the best performing model.
Helper functions
To help score the models, I have written a couple of helper functions. A takes a fitted model, predictions and the desired horizon, and it generates plots of actual, predicted and confidence intervals and the ASE for the given model.
eval_model <- function(response, predictions, pred_ul = NA, pred_ll = NA, model_name, AIC_val = 0, ending_point = length(response)) {
num_predictions = length(predictions)
test_stop <- length(response)
test_start <- test_stop - num_predictions + 1
compare_stop <- test_start - 1
compare_start <- compare_stop - num_predictions + 1
ASE <- mean((predictions - response[test_start:test_stop])^2)
# Build predictions dataframe
df <- data.frame('Predicted' = predictions)
df$Day = row(df)
df$Actual = response[test_start:test_stop]
df <- gather(df, key='Type', value, -Day)
#if we have enough data to plot num_predictions * 2, do it, else use num_predictions
starting_point <- compare_start - num_predictions + 1
plot_start <- ifelse(starting_point < 0, compare_start, starting_point)
day_multiplier <- ifelse(starting_point < 0, -1,-2)
# Build predicted vs actual dataframe
df <- rbind(df,
data.frame("Day"=c(((num_predictions-1)*day_multiplier):0),
"Type" = 'Actual',
value = response[plot_start:compare_stop]))
# Built UL/LL dataframes
ul <- data.frame("Day"=c(1:num_predictions), pred_ul)
ll <- data.frame("Day"=c(1:num_predictions), pred_ll)
# Build Plot
comparison_plot <- ggplot() +
geom_line(data=df, aes(Day + ending_point - num_predictions, value, color=Type)) +
geom_point(size=.75) +
labs(title=paste(model_name, 'Performance Evaluation'),
subtitle=paste0(num_predictions,'-Day Forecast'),
x='Day',
y='Positivity Rate',
caption=paste0('ASE: ',round(ASE,6),
'\nAIC: ',round(AIC_val,6)))
# Add confidence intervals if supplied
if (length(pred_ul) == length(predictions)){
comparison_plot = comparison_plot +
geom_line(aes(ul$Day + ending_point - num_predictions, ul$pred_ul),
color='grey70', linetype = "dashed")
}
if (length(pred_ll) == length(predictions)){
comparison_plot = comparison_plot +
geom_line(aes(ll$Day + ending_point - num_predictions, ll$pred_ll),
color='grey70', linetype = "dashed")
}
return(comparison_plot)
}########## example: ARIMA(12,2) ##########
e <- est.arma.wge(us$positivity, p=12, q=2)
preds <- fore.aruma.wge(us$positivity, phi = e$phi, theta=e$theta, n.ahead = 7, lastn = T, limits = F)
eval_model(us$positivity,preds$f, preds$ul, preds$ll, 'ARMA(12,2)', AIC_val = e$aic)
preds <- fore.aruma.wge(us$positivity, phi = e$phi, theta=e$theta, n.ahead = 90, lastn = T, limits = F)
eval_model(us$positivity,preds$f, preds$ul, preds$ll, 'ARMA(12,2)', AIC_val = e$aic) Another function calculates the rolling ASE, which is a measure of the ASE across several windows in time, for a fitted model. This function uses the evaluation function to plot each window and capture score, and then stores them in a list. This list is then averaged to find the average ASE across several time periods. By default I use a month of training data to evaluate a seven day period. The function finally returns a plot of all windows and actual values to visualize performance across the range of data.
rolling_ASE <- function (df, fitted_model, d=0, s=0, horizon=7, training_size=30, model_name, model_type = 'ARUMA', p, df_XDF=NA){
ASE = list(ASE = c(), plots = c(), multiplot = NA)
comp_df <- df %>% dplyr::select(date, positivity)
comp_df$preds = NA
names(comp_df) = c('date','Actual','Predicted')
test_stop <- length(df$positivity)
loop_end <- floor(test_stop/(training_size+horizon))
for (x in 1:loop_end){
test_start <- test_stop - horizon + 1
train_start <- test_start - training_size
train_stop <- test_start - 1
print(paste0('test window: ',test_start,':',test_stop,
', train window: ',train_start,':',train_stop))
data_window <- df$positivity[train_start:train_stop]
if(model_type == 'ARUMA') {
preds <- fore.aruma.wge(data_window,
phi=fitted_model$phi,
theta=fitted_model$theta,
s=s,
d=d,
n.ahead = horizon,
lastn = F,
limits = F)
pred_object <- preds$f
}
if(model_type == 'SigPlusNoise') {
preds <- fore.sigplusnoise.wge(data_window, max.p = p, n.ahead = horizon, limits=F)
pred_object <- preds$f
}
if(model_type == 'NNFOR') {
ts_la <- ts(data_window, start = '1')
mlp_model = mlp(ts_la, lags = horizon, hd.auto.type = 'cv')
?mlp
preds <- predict(mlp_model, horizon)
preds$ul <- NA
preds$ll <- NA
pred_object <- as.numeric(preds$mean)
}
if(model_type == 'VAR') {
vfit=VAR(cbind(Positivity = df$positivity, df_XDF)[train_start:train_stop,], p=p, type='both', season = s)
preds=predict(vfit,n.ahead=7)
pred_object <- preds$fcst$Positivity[,1]
preds$ul <- preds$fcst$Positivity[,3]
preds$ll <- preds$fcst$Positivity[,2]
}
a <- mean((pred_object - df$positivity[test_start:test_stop])^2)
print(paste('Window ASE:', a))
ASE$ASE[x] <- a
comp_df$Predicted[test_start:test_stop] = pred_object
ASE$plots[x] <-
plot(eval_model(
data_window,
pred_object,
model_name = model_name,
AIC_val = ifelse(model_type == 'ARUMA', fitted_model$aic, 0),
pred_ul = preds$ul,
pred_ll = preds$ll,
ending_point = test_stop))
test_stop = test_stop - training_size
}
ASE$multiplot <- plot(gather(comp_df, key = Type, value = 'Positivity_Rate', -date) %>%
ggplot(aes(date, Positivity_Rate, color=Type)) + geom_line() +
labs(
title=paste(model_name, 'Performance Evaluation'),
subtitle=paste0(horizon,'-Day Forecast Rolling Window'),
x='Day',
y='Positivity Rate',
caption=paste0('Mean ASE: ',round(mean(ASE$ASE),6))
)
)
return(ASE)
}########## example: ARIMA(12,2) ##########
e <- est.arma.wge(us$positivity, p=12, q=2)
test <- rolling_ASE(us, e, s=12, d=1, horizon=7, model_name = 'ARMA(12,2)')Univariate Analysis
The steps for building and evaluating our models will proceed as follows: - Plot and examine the data - Realization - Spectral Density - Autocorrelation Function - Evaluate stationarity assumptions. There are three stationary conditions: - Mean does not depend on time - Variance is finite and does not depend on time - Correlation only depends on how far apart observations are in time, not where they are - Build models - Check that serial correlation has been removed (AIC scores, ACF, PACF plots) - Check residuals for white noise - Evaluate model forecast performance
Louisiana
Plot the Data:
# Plot the data
plotts.sample.wge(la$positivity)
## $autplt
## [1] 1.0000000 0.9393071 0.8796384 0.8605446 0.8399214 0.8267488 0.8490105
## [8] 0.8610586 0.8121303 0.7578559 0.7277727 0.6986874 0.6752319 0.6832333
## [15] 0.6821986 0.6292681 0.5668868 0.5319428 0.4966430 0.4664633 0.4593505
## [22] 0.4569352 0.4031459 0.3462679 0.3230427 0.2910256
##
## $freq
## [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
## [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
## [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
## [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
## [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
## [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
## [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
## [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
## [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
## [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
## [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
## [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
## [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
## [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
## [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
## [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
## [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
##
## $db
## [1] 13.9602603 17.5611464 11.9599796 7.6668080 6.4697337 1.6128896
## [7] -2.8203865 -5.0586421 -9.6952719 -7.3153505 -23.1165770 -8.5853171
## [13] -6.4269213 -10.4512454 -3.9370940 -6.6895400 -5.0361316 -8.8174379
## [19] -8.8408129 -7.7524816 -12.7800934 -10.8831615 -12.3443173 -7.4187515
## [25] -13.2201935 -12.5483706 -12.6072411 -10.4639780 -10.3061541 -5.7896562
## [31] -6.9527205 -4.5729994 -2.6565916 -11.1992997 4.9222172 -3.5726203
## [37] -17.5737307 -6.4938454 -6.9104502 -10.1859592 -7.9690657 -14.6326883
## [43] -12.8126497 -12.2660625 -9.3597016 -18.1927808 -10.9197193 -15.3687381
## [49] -12.2455703 -14.3938321 -14.3987676 -8.0215954 -11.2276907 -13.8151299
## [55] -18.4209672 -14.1462588 -13.9178056 -10.3542921 -18.1090269 -9.4175468
## [61] -11.6144991 -12.6302015 -16.0793212 -17.8525413 -14.4047311 -13.6480654
## [67] -19.4031903 -6.2595560 -6.4716425 -0.1702762 -7.4453558 -16.3967146
## [73] -11.3842620 -12.6276061 -23.1800642 -15.0448944 -12.1965307 -17.0067301
## [79] -13.2831602 -15.2489091 -11.9935760 -10.7178034 -13.0921498 -25.3789882
## [85] -14.9455352 -23.1102712 -24.1520012 -20.1648753 -14.8792309 -16.2779671
## [91] -15.7492788 -11.4820211 -17.6783892 -34.5124466 -17.2522571 -14.1063598
## [97] -22.9668371 -28.1216282 -22.5894863 -17.0841672 -18.0459368 -11.7319206
## [103] -14.1578096 -16.1153463 -14.0821580 -12.1054220 -21.5900996 -29.1884372
## [109] -31.9989494 -22.2731413 -16.8603667 -42.1727225 -30.0127617 -10.3486417
## [115] -15.6342856 -16.8323448 -23.1519542 -18.4824601 -20.6476432 -22.7906872
## [121] -16.9156395 -31.3160094 -23.1540488
##
## $dbz
## [1] 12.566881 12.293334 11.835242 11.189378 10.351358 9.315923
## [7] 8.077572 6.631985 4.979074 3.129376 1.116753 -0.979425
## [13] -3.007511 -4.747472 -6.005907 -6.756174 -7.136643 -7.314196
## [19] -7.399551 -7.449649 -7.491625 -7.533833 -7.564111 -7.546353
## [25] -7.426789 -7.156016 -6.718379 -6.146017 -5.505068 -4.868567
## [31] -4.297108 -3.833194 -3.503840 -3.325433 -3.307926 -3.457655
## [37] -3.778856 -4.274127 -4.943857 -5.784457 -6.784834 -7.920367
## [43] -9.144086 -10.377572 -11.511317 -12.430793 -13.068592 -13.440575
## [49] -13.624521 -13.706896 -13.748785 -13.782319 -13.821469 -13.872769
## [55] -13.939776 -14.020358 -14.098741 -14.136324 -14.068719 -13.820057
## [61] -13.338769 -12.633433 -11.773842 -10.856577 -9.969285 -9.175272
## [67] -8.514349 -8.009317 -7.672104 -7.507927 -7.517585 -7.698307
## [73] -8.043468 -8.541337 -9.172968 -9.909643 -10.711031 -11.526432
## [79] -12.302013 -12.994333 -13.583945 -14.078851 -14.503984 -14.884222
## [85] -15.231765 -15.543263 -15.806171 -16.010199 -16.157419 -16.264982
## [91] -16.359111 -16.464519 -16.594649 -16.745743 -16.895822 -17.009646
## [97] -17.050528 -16.996280 -16.851085 -16.645161 -16.423033 -16.229091
## [103] -16.097778 -16.049668 -16.091173 -16.215626 -16.404893 -16.631990
## [109] -16.865795 -17.078322 -17.253005 -17.390290 -17.507267 -17.631541
## [115] -17.792797 -18.015496 -18.313830 -18.688042 -19.120778 -19.573264
## [121] -19.984024 -20.276769 -20.383632# Look deeper at spectral density
parzen.wge(la$positivity, trun=100)## $freq
## [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
## [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
## [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
## [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
## [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
## [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
## [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
## [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
## [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
## [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
## [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
## [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
## [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
## [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
## [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
## [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
## [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
##
## $pzgram
## [1] 14.8637790 14.2011980 12.6485602 10.0671892 6.8101322 3.6227123
## [7] 0.7036374 -2.0633900 -4.3842228 -6.5433414 -8.2611328 -8.4389751
## [13] -7.6413777 -6.8009486 -6.2911189 -6.1602063 -6.3694419 -6.8987848
## [19] -7.7062083 -8.6863538 -9.5077432 -9.6557081 -9.2342817 -8.9830525
## [25] -9.2278061 -9.6798562 -9.7913461 -9.1393431 -7.8987016 -6.5685260
## [31] -5.5392348 -4.6515485 -2.9494978 -0.8164328 0.2786933 -0.3491239
## [37] -2.7974368 -6.3739983 -8.8812547 -10.0296933 -11.2356596 -12.5004769
## [43] -13.5470125 -13.8931208 -13.8616074 -14.0405632 -14.3322121 -14.6024260
## [49] -14.6257072 -13.9303719 -12.8233618 -12.2324784 -12.6920818 -14.1670234
## [55] -15.7658464 -15.9690501 -15.1426719 -14.4520904 -14.0367957 -13.7506921
## [61] -13.7890897 -14.4329786 -15.5991766 -16.6962898 -16.3502735 -14.2029150
## [67] -11.2635191 -7.8543781 -5.0149509 -3.6453479 -4.0326645 -6.1283961
## [73] -9.2993396 -11.9398547 -13.3944641 -13.8986525 -13.5298106 -13.2516942
## [79] -13.2951733 -13.1402295 -12.7646338 -12.8686582 -13.8904306 -15.5765112
## [85] -17.3015466 -18.7823532 -19.4344998 -18.5676263 -16.9718976 -15.3765724
## [91] -14.3372998 -14.4082388 -15.7083639 -17.3837409 -18.0321588 -18.3635612
## [97] -19.4736057 -20.5164707 -19.8122370 -18.0227462 -16.4113405 -15.4494854
## [103] -14.9042987 -14.1846154 -13.5910594 -14.0910964 -16.3214414 -20.2716849
## [109] -23.0386491 -22.0845605 -20.6856531 -18.5019235 -16.0286579 -14.7041570
## [115] -15.0159308 -16.8615493 -19.4498201 -21.0946956 -21.2866714 -21.0455908
## [121] -21.2006319 -22.0250180 -22.6403932
There is strong evidence of serial correlation in this data as can be seen in the peaks in the spectral density plots and the slowly damping, sinusoidal behavior in the autocorrelation function plots. There is significant evidence of a weekly trend, as can be seen in the Parzen Window (spectral density plot). The frequency peak at .14 indicates a weekly trend (1/7 = .143), with potentially an echo at .28. Increasing the truncation point of the Parzen function helps highlight some additional points that suggest possibly a monthly factor as well.
Evaluate Stationarity:
There are several indications that suggest the data is NOT stationary, the most obvious being the strong weekly trend with strong peaks each weekend, which effectively increases the mean on weekends. There is a difference in variance which can be seen in the early weeks of the pandemic that have tapered off as time progresses.
Model Construction:
I start by differencing the data, which filters long-term trending behavior and leaves cyclical patterns. Differencing once leaves a good bit of cyclical pattern in the ACF, especially at lag 7.
Identify correlation structures:
Here we are going to look at different filters to try to model out some of the correlation in the data. We’ll also look at a factor table to see if we can identify a weekly or monthly trend.
The overfit table shows strong evidence of 1-B, 1+.445B+B^2, 1-1.247B+B^2 and 1+.802B+B^2 being present in the model. This is more than enough evidence to proceed with the weekly seasonal model.
# Take out 1-B
la.d1 <- artrans.wge(la$positivity,1)
# Taking out another 1-B.. looks like it isn't enough
artrans.wge(la.d1,1)
## [1] -0.12001810013 0.11197961593 0.01148788128 -0.09768340278 0.08729785421
## [6] -0.01424222648 -0.00245022335 0.02856748199 0.03663734739 -0.03020620067
## [11] -0.08490402636 0.07586325598 -0.04411939164 0.01405163464 0.00603443239
## [16] 0.05326455229 -0.07040279424 -0.01703609861 -0.00895351037 -0.00138728940
## [21] 0.03352680899 -0.02059646162 0.08058041871 -0.01655174992 -0.12569990013
## [26] 0.07715699498 0.00970828146 -0.02828671007 0.04922032309 -0.00137427599
## [31] 0.01757626942 -0.08886926099 0.04368065913 -0.01926462275 0.01546202308
## [36] 0.01132636604 0.01489236125 0.02424383991 -0.09006170447 0.04445552513
## [41] 0.01603133393 -0.01939040709 0.00042712593 0.01826882185 -0.00377191437
## [46] -0.01780837066 -0.00002107674 0.02992415119 -0.03931100761 0.02461497660
## [51] 0.01422652210 -0.01024811730 -0.03048764593 0.01974597812 0.00900167078
## [56] 0.00819385113 -0.03064462453 0.00950030564 0.04636873235 -0.05359735437
## [61] 0.00719330921 0.00019993853 0.02416632150 -0.01474365246 -0.01262301520
## [66] 0.04615993257 -0.05738147393 0.00585647802 0.02531654811 -0.00093831949
## [71] -0.01817074820 0.04329415542 -0.04213960697 0.00626173920 -0.00566945709
## [76] 0.00745556817 0.00195300117 0.00524211825 0.00262408214 0.00460223860
## [81] -0.02713178495 0.00522217386 0.01741746191 -0.01115368628 0.00752134400
## [86] 0.00774049302 0.01478958929 -0.06086352002 0.03782370394 -0.00571849966
## [91] 0.00319051199 0.00398135599 0.04963294085 -0.05888886023 -0.03562008158
## [96] 0.04411647176 0.00224009636 0.00153450585 -0.00789657627 0.03008822968
## [101] -0.01229604121 -0.05886386305 0.03469231074 0.00710503644 0.00257321289
## [106] 0.00022733426 0.04154065525 -0.04786083785 -0.01709010640 0.01517072373
## [111] 0.00760007924 0.00649024584 0.02373510913 -0.00284543631 -0.01272601102
## [116] -0.08239005417 0.04713628416 0.02718539078 -0.00738986440 -0.00059194356
## [121] 0.03070118614 -0.02070185623 -0.05368870817 0.02827825825 0.01409853916
## [126] 0.00576960680 -0.00084032111 0.02024613738 -0.01441860235 -0.03524539380
## [131] 0.01209539642 0.01070138410 0.00084646748 0.00828549360 0.01419305283
## [136] -0.00792484179 -0.04416309640 0.00375735815 0.03480591454 -0.00913183129
## [141] 0.00713310866 0.01075547427 -0.01421246068 -0.03455009149 0.02639598188
## [146] 0.00021892746 0.01219897879 -0.01321799746 0.03121307054 -0.01499341370
## [151] -0.04776882472 0.02695663191 0.00674473559 0.01212144857 -0.00504886828
## [156] 0.03724066022 -0.06437056017 -0.00005126073 0.02327608778 0.00336498625
## [161] -0.00103675954 0.00310067789 0.01971938490 -0.02913216776 -0.01216822883
## [166] 0.01335914467 -0.00013283852 0.00748879028 0.01075283531 -0.01444866773
## [171] 0.01637641523 -0.04955039165 0.02621505529 0.01434378900 -0.01189553928
## [176] 0.00512143710 0.01253283708 -0.00618458919 -0.03496356980 0.03070706559
## [181] -0.00885731849 0.00400893776 0.00250957556 0.03818324918 -0.03313117378
## [186] -0.04599125522 0.04620206871 -0.00982100992 0.00260166582 0.00312275101
## [191] 0.02014836905 -0.02033125343 -0.02344872791 0.01835402238 0.01194806885
## [196] -0.00898995412 0.00283657842 0.02190296124 -0.02927985203 -0.01345968115
## [201] 0.01367358904 0.00621727754 0.00114851488 0.00516535494 0.01112481513
## [206] -0.01896296615 -0.02146644421 0.01711274737 0.01075193799 -0.00423858946
## [211] -0.00013615570 0.04181976478 -0.04444093558 -0.03623206054 0.03222659037
## [216] 0.01219654120 -0.00589440433 0.00152342606 0.02251320599 -0.02345323556
## [221] -0.03253380372 0.02430159488 0.01522738477 -0.00666540167 0.00040987785
## [226] 0.01831575860 -0.01954011115 -0.02892179674 0.03014970836 0.00318487094
## [231] 0.00191084383 -0.00973136741 0.03725493099 -0.03444842445 -0.04309372782
## [236] 0.04491734584 0.01680699336 -0.01430450447 -0.00035522384 0.02947585241
## [241] -0.00995478282 -0.07682562404 0.04821360473 0.01244096381# strong evidence of a 7 day cycle + possibly a monthly cycle
la.s7 <- artrans.wge(la$positivity, c(rep(0,6),1))
# Add (1-B) + s7
la.d1.s7 <- artrans.wge(la.s7,1)
# Try to model weekly + monthly data
la.s7.s12 <- artrans.wge(la$positivity, c(rep(0,6),1,0,0,0,0,1))
# Overfit table
est.ar.wge(la$positivity, p=15)##
## Coefficients of Original polynomial:
## 0.7229 0.1790 0.2398 -0.1301 -0.0733 0.0620 0.3964 -0.3014 -0.2877 -0.1168 0.0513 0.0441 0.0674 0.3147 -0.1753
##
## Factor Roots Abs Recip System Freq
## 1-0.9959B 1.0041 0.9959 0.0000
## 1-1.2348B+0.9762B^2 0.6324+-0.7902i 0.9881 0.1426
## 1+0.4517B+0.9649B^2 -0.2340+-0.9908i 0.9823 0.2869
## 1-1.8259B+0.8957B^2 1.0193+-0.2785i 0.9464 0.0424
## 1+1.6899B+0.8887B^2 -0.9508+-0.4704i 0.9427 0.4269
## 1+1.1264B+0.7953B^2 -0.7082+-0.8694i 0.8918 0.3588
## 1+0.8673B -1.1530 0.8673 0.5000
## 1-0.2981B+0.6757B^2 0.2206+-1.1964i 0.8220 0.2210
## 1-0.5036B 1.9857 0.5036 0.0000
##
##
##
## Coefficients of Original polynomial:
## 0.7229 0.1790 0.2398 -0.1301 -0.0733 0.0620 0.3964 -0.3014 -0.2877 -0.1168 0.0513 0.0441 0.0674 0.3147 -0.1753
##
## Factor Roots Abs Recip System Freq
## 1-0.9959B 1.0041 0.9959 0.0000
## 1-1.2348B+0.9762B^2 0.6324+-0.7902i 0.9881 0.1426
## 1+0.4517B+0.9649B^2 -0.2340+-0.9908i 0.9823 0.2869
## 1-1.8259B+0.8957B^2 1.0193+-0.2785i 0.9464 0.0424
## 1+1.6899B+0.8887B^2 -0.9508+-0.4704i 0.9427 0.4269
## 1+1.1264B+0.7953B^2 -0.7082+-0.8694i 0.8918 0.3588
## 1+0.8673B -1.1530 0.8673 0.5000
## 1-0.2981B+0.6757B^2 0.2206+-1.1964i 0.8220 0.2210
## 1-0.5036B 1.9857 0.5036 0.0000
##
## ## $phi
## [1] 0.72293214 0.17904609 0.23984475 -0.13008877 -0.07327756 0.06198856
## [7] 0.39641995 -0.30136793 -0.28766989 -0.11683003 0.05126236 0.04414203
## [13] 0.06741100 0.31465706 -0.17526686
##
## $res
## [1] -0.03093945037 0.02316236956 -0.05435390974 -0.00502183438 0.00893904094
## [6] -0.01985489760 0.02518009493 -0.00323618159 -0.01080959161 0.05450484314
## [11] 0.05434986643 0.01894490423 -0.02611320488 0.00463390908 -0.02988421379
## [16] -0.01595229496 0.00519612379 0.02334120053 -0.01637060949 0.02673726934
## [21] -0.02746859385 -0.02094560817 0.00948456862 -0.02138901615 0.01445368950
## [26] 0.03925759433 -0.03524684298 -0.01054355581 0.00942017257 -0.01413700563
## [31] 0.00694955701 -0.04274641357 0.00074863942 -0.01905898380 0.00565668039
## [36] -0.02889700609 -0.01189626903 -0.01409257323 -0.02753259822 0.00034614558
## [41] -0.00850118185 0.00066293073 0.01340145691 0.00903591305 -0.02212165781
## [46] -0.02874659746 -0.03762965113 0.00770188222 -0.00551425507 0.01429056095
## [51] -0.01787291268 -0.00059923134 -0.00156357314 -0.00991733273 -0.01161547248
## [56] -0.00896680848 -0.01191706546 0.01956372890 -0.01415446599 -0.02954068870
## [61] 0.01749759469 0.00517144936 0.00130098304 -0.02678301296 0.00678500035
## [66] 0.00814910002 -0.01142255741 0.00807402729 -0.01465019447 -0.01094130636
## [71] 0.00186169458 -0.00131120746 -0.00604921370 0.03148190920 -0.02682804486
## [76] 0.00412801364 -0.00174523049 -0.00123813447 -0.01123950065 0.00317430085
## [81] -0.00282394477 0.00263284700 -0.00311323827 -0.00067644099 0.00368140353
## [86] -0.00288795943 0.00256798041 -0.00542778472 0.02242509916 -0.02085880014
## [91] 0.00566567188 -0.00770799376 0.00408359714 0.00116957504 0.04409157264
## [96] -0.00638623605 -0.02161682831 -0.00287760560 0.00660384509 0.01811151176
## [101] -0.00261361664 -0.00116705444 0.01225968475 0.00024633298 -0.00248237358
## [106] -0.00568342337 0.00470324489 0.00321782831 0.01666260609 -0.00190960440
## [111] 0.01333832066 0.00304272435 0.00417870431 0.00846533013 0.02955830254
## [116] 0.00731876833 0.01483670557 -0.03764748296 -0.01929639747 0.00351744585
## [121] 0.01039928781 -0.00854192500 0.00054005230 0.01490819192 0.00938605645
## [126] 0.00770159066 -0.00417226682 0.00046251487 -0.00964401877 -0.00644861392
## [131] -0.00429676443 0.01314316916 0.00481705246 -0.00584066542 -0.00678062295
## [136] -0.00184183800 -0.00395409583 0.00431774407 -0.00443480516 -0.01568071375
## [141] 0.00408165357 0.00177380610 0.00267825774 -0.00603350603 -0.00899924471
## [146] -0.01134564320 0.00916042653 -0.00580768182 0.00418374479 -0.01636468317
## [151] 0.00200478699 0.00222983810 -0.01240003293 -0.00612019352 -0.01237420014
## [156] 0.00393755026 -0.00059175998 0.01909125446 -0.02822531155 -0.00629549502
## [161] 0.00136862439 0.00741735403 -0.00360759578 -0.00595179654 -0.00683535968
## [166] -0.00068930161 0.00915453494 0.00076502854 -0.00714786447 -0.00200711739
## [171] 0.00654471037 -0.02202500232 0.01806906077 -0.01393317717 -0.00412223383
## [176] 0.00358651115 -0.00263152438 -0.00939463150 -0.00261176383 0.00370279100
## [181] -0.00694669263 0.00588280941 -0.00983992258 -0.00498181553 -0.00977576459
## [186] 0.02602126906 0.00301226969 -0.01863472947 -0.00355022985 -0.00860476524
## [191] 0.00329891257 -0.00487912168 -0.00550424536 -0.00582125658 0.00506704835
## [196] -0.00146964615 0.00898300503 0.00374046125 -0.00121913464 -0.00277754085
## [201] -0.00943344986 0.00076418798 -0.00672419290 -0.00331745977 0.00225067615
## [206] 0.00446921729 0.00164313697 0.00131682864 -0.00458535259 -0.00253809423
## [211] 0.00121465831 0.00211380925 -0.00523795138 0.02301279532 0.00714116899
## [216] -0.01666709162 -0.00992813954 0.00061765639 0.00485997939 -0.00195597126
## [221] -0.00276673410 0.00177202064 -0.00053047284 0.00152742254 0.00312650034
## [226] 0.00371745232 -0.00117132795 -0.00724414542 -0.00125482227 -0.00101012899
## [231] 0.00733384230 -0.00006790254 0.00697149684 -0.00651165092 0.01388212638
## [236] 0.00286561357 -0.01701313400 -0.00316294895 0.01119596746 0.01067369287
## [241] 0.00460069718 0.01018875659 0.02827800653 -0.01218119477 -0.00701711434
## [246] -0.01030342246
##
## $avar
## [1] 0.0002103418
##
## $aic
## [1] -8.336695
##
## $aicc
## [1] -7.317654
##
## $bic
## [1] -8.108706# Generate a factor table for a (1-B^7)
tswge::factor.wge(c(rep(0,6),1))##
## Coefficients of Original polynomial:
## 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000
##
## Factor Roots Abs Recip System Freq
## 1-1.0000B 1.0000 1.0000 0.0000
## 1+0.4450B+1.0000B^2 -0.2225+-0.9749i 1.0000 0.2857
## 1-1.2470B+1.0000B^2 0.6235+-0.7818i 1.0000 0.1429
## 1+1.8019B+1.0000B^2 -0.9010+-0.4339i 1.0000 0.4286
##
## Estimate parameters and check residuals:
TSWGE offers an aic5 function that iterates through values of P (phi) and Q (theta), fits a model, and returns the top 5 performing models by either AIC or BIC value. The AIC and BIC penalize models for the number of parameters used, with BIC favoring simpler models even more so than AIC. This gives us a pretty good idea of the best fitting model as a starting point.
The ACF and PACF plots look like white noise, as expected. The Ljung-Box test fails to reject the null hyphothesis that the residuals are white noise. The generated spectral densities look like they have a fairly good fit, the final peak is overpronounced in our generations compared to the actual data. The ACF comparisons are all over the board, but with peaks in the same general spots (lags at 7, 14, 21). The generated realizations look quite a bit different, but the spectral densities and ACF plots look fairly similar, indicating this could be a useful model.
# Estimate the model params
aic5.wge(la.s7, type = 'bic', p=0:20,q=0:2) #15## ---------WORKING... PLEASE WAIT...
##
##
## Error in aic calculation at 18 1
## Five Smallest Values of bic## p q bic
## 46 15 0 -8.256251
## 49 16 0 -8.249310
## 48 15 2 -8.247515
## 42 13 2 -8.246550
## 52 17 0 -8.242890e <- est.arma.wge(la.s7, p=15, q=0)##
## Coefficients of Original polynomial:
## 0.6819 0.2329 0.1795 -0.1438 0.0353 0.0725 -0.6438 0.3639 0.0825 -0.0664 -0.1285 0.0148 0.2163 -0.3415 0.3063
##
## Factor Roots Abs Recip System Freq
## 1-1.8995B+0.9868B^2 0.9625+-0.2950i 0.9934 0.0473
## 1+0.9661B+0.9242B^2 -0.5226+-0.8994i 0.9614 0.3338
## 1+1.3916B+0.9115B^2 -0.7634+-0.7172i 0.9547 0.3800
## 1+1.8666B+0.8970B^2 -1.0405+-0.1795i 0.9471 0.4728
## 1-0.9455B 1.0576 0.9455 0.0000
## 1-0.1935B+0.8186B^2 0.1182+-1.0989i 0.9048 0.2329
## 1-1.3104B+0.7391B^2 0.8865+-0.7531i 0.8597 0.1121
## 1-0.5571B+0.7181B^2 0.3879+-1.1145i 0.8474 0.1967
##
## # Ljung Box Test shows white noise residuals
ljung.wge(artrans.wge(la.s7, phi.tr = e$phi))## Obs 0.03251153 0.005953455 -0.04811914 0.1123665 0.05091219 -0.006740414 0.02838045 0.05305956 -0.01653444 0.03643441 0.01285373 0.02944821 0.008132136 -0.003795011 0.08950678 -0.02009765 0.05161991 0.02098249 0.0549419 0.02685266 -0.081336 0.07919188 -0.08543065 0.0593491## $test
## [1] "Ljung-Box test"
##
## $K
## [1] 24
##
## $chi.square
## [1] 15.45967
##
## $df
## [1] 24
##
## $pval
## [1] 0.9065855ljung.wge(artrans.wge(la.s7, phi.tr = e$phi), K = 48)
## Obs 0.03251153 0.005953455 -0.04811914 0.1123665 0.05091219 -0.006740414 0.02838045 0.05305956 -0.01653444 0.03643441 0.01285373 0.02944821 0.008132136 -0.003795011 0.08950678 -0.02009765 0.05161991 0.02098249 0.0549419 0.02685266 -0.081336 0.07919188 -0.08543065 0.0593491 -0.03150929 0.05970946 -0.09874056 0.04136814 0.06256806 -0.0520347 -0.04590716 0.06908596 0.04504184 -0.05018569 -0.09408055 -0.05452445 0.1304085 -0.04294655 -0.03968746 -0.02780619 0.0549043 -0.141726 0.1096732 -0.03633432 -0.02496781 -0.03631315 0.05437685 -0.02354808## $test
## [1] "Ljung-Box test"
##
## $K
## [1] 48
##
## $chi.square
## [1] 44.94805
##
## $df
## [1] 48
##
## $pval
## [1] 0.5986569acf(e$res)
pacf(e$res)
dev.off()## null device
## 1#Compare Spectral Densities
sims = 5
SpecDen = parzen.wge(la$positivity, plot = "FALSE")
plot(SpecDen$freq, SpecDen$pzgram, type = "l", lwd = 6)
for( i in 1: sims)
{
SpecDen2 = parzen.wge(gen.aruma.wge(246,s = 7, phi = e$phi, plot ="FALSE"), plot = "FALSE")
lines(SpecDen2$freq,SpecDen2$pzgram, lwd = 2, col = "red")
}
#Compare ACFs
sims = 5
ACF = acf(la$positivity, plot = "FALSE")
plot(ACF$lag ,ACF$acf , type = "l", lwd = 6)
for( i in 1: sims)
{
ACF2 = acf(gen.aruma.wge(246, s = 7, phi = e$phi, plot = "FALSE"), plot = "FALSE")
lines(ACF2$lag ,ACF2$acf, lwd = 2, col = "red")
}
#Compare Generated Realizations
eGen = gen.aruma.wge(246, s = 7, phi = e$phi, vara = e$avar)
plotts.sample.wge(eGen)## $autplt
## [1] 1.0000000 0.9736991 0.9626871 0.9439085 0.9248900 0.9046038 0.8813618
## [8] 0.8742908 0.8401772 0.8238030 0.8047006 0.7929718 0.7858042 0.7776364
## [15] 0.7873074 0.7724358 0.7735335 0.7716518 0.7714489 0.7694762 0.7618503
## [22] 0.7669221 0.7428109 0.7309969 0.7125539 0.6936437
##
## $freq
## [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
## [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
## [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
## [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
## [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
## [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
## [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
## [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
## [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
## [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
## [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
## [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
## [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
## [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
## [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
## [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
## [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
##
## $db
## [1] 18.959905 11.822348 8.431301 5.279660 2.998245 1.920224
## [7] 2.019772 0.889596 -4.071676 2.360229 3.035087 5.596055
## [13] -14.718606 -5.706822 -5.895885 -9.281965 -10.965889 -8.360391
## [19] -11.536625 -8.941074 -11.201632 -9.397568 -10.128222 -11.581647
## [25] -10.530920 -10.415388 -12.223071 -13.714795 -15.929249 -11.416607
## [31] -14.650649 -13.000973 -8.652443 -6.858406 -9.184427 -13.538787
## [37] -13.512267 -13.838175 -11.400413 -13.597119 -16.074031 -17.743589
## [43] -19.187122 -15.621058 -16.246249 -19.156378 -20.689526 -21.787558
## [49] -21.207093 -23.609956 -17.776729 -19.532000 -15.296678 -16.045822
## [55] -15.673240 -20.973239 -18.562394 -17.623196 -16.395773 -16.102934
## [61] -13.468320 -18.070707 -17.961197 -22.573351 -15.967344 -17.201750
## [67] -16.657263 -11.192350 -9.885698 -14.675763 -8.173531 -11.008147
## [73] -18.271583 -31.978990 -23.808425 -24.623456 -29.502731 -24.997538
## [79] -17.720580 -23.963103 -14.550412 -15.840960 -17.610135 -26.105177
## [85] -24.648677 -19.517523 -17.956652 -21.135251 -25.198270 -20.907685
## [91] -23.646234 -19.497703 -25.984641 -20.813359 -28.070662 -21.818354
## [97] -29.093402 -18.943961 -33.424419 -24.213975 -21.867636 -25.267836
## [103] -20.129915 -18.508756 -3.395199 -15.594638 -9.780505 -12.382323
## [109] -23.644788 -23.055684 -20.305960 -17.995426 -22.721691 -22.945210
## [115] -25.848928 -15.746649 -15.911528 -19.473489 -26.933756 -22.765069
## [121] -17.658643 -23.506246 -22.314047
##
## $dbz
## [1] 12.98338914 12.67366873 12.15552791 11.42660225 10.48495692
## [6] 9.33124119 7.97286582 6.43161294 4.75581155 3.03426865
## [11] 1.39679121 -0.02681045 -1.18680272 -2.15349125 -3.06123556
## [16] -4.02262187 -5.08250916 -6.20823593 -7.30055818 -8.23754381
## [21] -8.95433325 -9.49493285 -9.97289978 -10.48822020 -11.07456247
## [26] -11.68590510 -12.21180282 -12.53019870 -12.58691479 -12.43198324
## [31] -12.17306092 -11.90820127 -11.69804179 -11.57226267 -11.54552778
## [36] -11.62965214 -11.83896601 -12.19009842 -12.69849368 -13.37343017
## [41] -14.21210568 -15.19254148 -16.26576559 -17.35119214 -18.34460413
## [46] -19.14695255 -19.70176485 -20.00808893 -20.09903760 -20.01672232
## [51] -19.80703194 -19.52417350 -19.22771230 -18.97079596 -18.78820767
## [56] -18.68937058 -18.65576573 -18.64168561 -18.58055500 -18.40142305
## [61] -18.05518080 -17.53712130 -16.88829227 -16.17555161 -15.46788247
## [66] -14.82260989 -14.28210319 -13.87578779 -13.62336632 -13.53747101
## [71] -13.62529072 -13.88913125 -14.32585590 -14.92501128 -15.66537773
## [76] -16.51005025 -17.40168428 -18.26289651 -19.00989574 -19.58185503
## [81] -19.96951744 -20.21667709 -20.39246380 -20.55966131 -20.75796337
## [86] -21.00171502 -21.28492705 -21.58899094 -21.89083122 -22.16819344
## [91] -22.39672151 -22.53526310 -22.50600964 -22.19544773 -21.50929252
## [96] -20.45624686 -19.16111365 -17.78777125 -16.46736359 -15.28023881
## [101] -14.26794488 -13.44900781 -12.82996291 -12.41168250 -12.19267076
## [106] -12.17060484 -12.34282141 -12.70602136 -13.25515670 -13.98119877
## [111] -14.86726515 -15.88261404 -16.97503703 -18.06562416 -19.05583463
## [116] -19.85709990 -20.43230235 -20.81111603 -21.06017500 -21.23853777
## [121] -21.37384042 -21.46458579 -21.49722005plotts.sample.wge(la$positivity)## $autplt
## [1] 1.0000000 0.9393071 0.8796384 0.8605446 0.8399214 0.8267488 0.8490105
## [8] 0.8610586 0.8121303 0.7578559 0.7277727 0.6986874 0.6752319 0.6832333
## [15] 0.6821986 0.6292681 0.5668868 0.5319428 0.4966430 0.4664633 0.4593505
## [22] 0.4569352 0.4031459 0.3462679 0.3230427 0.2910256
##
## $freq
## [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
## [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
## [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
## [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
## [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
## [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
## [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
## [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
## [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
## [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
## [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
## [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
## [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
## [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
## [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
## [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
## [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
##
## $db
## [1] 13.9602603 17.5611464 11.9599796 7.6668080 6.4697337 1.6128896
## [7] -2.8203865 -5.0586421 -9.6952719 -7.3153505 -23.1165770 -8.5853171
## [13] -6.4269213 -10.4512454 -3.9370940 -6.6895400 -5.0361316 -8.8174379
## [19] -8.8408129 -7.7524816 -12.7800934 -10.8831615 -12.3443173 -7.4187515
## [25] -13.2201935 -12.5483706 -12.6072411 -10.4639780 -10.3061541 -5.7896562
## [31] -6.9527205 -4.5729994 -2.6565916 -11.1992997 4.9222172 -3.5726203
## [37] -17.5737307 -6.4938454 -6.9104502 -10.1859592 -7.9690657 -14.6326883
## [43] -12.8126497 -12.2660625 -9.3597016 -18.1927808 -10.9197193 -15.3687381
## [49] -12.2455703 -14.3938321 -14.3987676 -8.0215954 -11.2276907 -13.8151299
## [55] -18.4209672 -14.1462588 -13.9178056 -10.3542921 -18.1090269 -9.4175468
## [61] -11.6144991 -12.6302015 -16.0793212 -17.8525413 -14.4047311 -13.6480654
## [67] -19.4031903 -6.2595560 -6.4716425 -0.1702762 -7.4453558 -16.3967146
## [73] -11.3842620 -12.6276061 -23.1800642 -15.0448944 -12.1965307 -17.0067301
## [79] -13.2831602 -15.2489091 -11.9935760 -10.7178034 -13.0921498 -25.3789882
## [85] -14.9455352 -23.1102712 -24.1520012 -20.1648753 -14.8792309 -16.2779671
## [91] -15.7492788 -11.4820211 -17.6783892 -34.5124466 -17.2522571 -14.1063598
## [97] -22.9668371 -28.1216282 -22.5894863 -17.0841672 -18.0459368 -11.7319206
## [103] -14.1578096 -16.1153463 -14.0821580 -12.1054220 -21.5900996 -29.1884372
## [109] -31.9989494 -22.2731413 -16.8603667 -42.1727225 -30.0127617 -10.3486417
## [115] -15.6342856 -16.8323448 -23.1519542 -18.4824601 -20.6476432 -22.7906872
## [121] -16.9156395 -31.3160094 -23.1540488
##
## $dbz
## [1] 12.566881 12.293334 11.835242 11.189378 10.351358 9.315923
## [7] 8.077572 6.631985 4.979074 3.129376 1.116753 -0.979425
## [13] -3.007511 -4.747472 -6.005907 -6.756174 -7.136643 -7.314196
## [19] -7.399551 -7.449649 -7.491625 -7.533833 -7.564111 -7.546353
## [25] -7.426789 -7.156016 -6.718379 -6.146017 -5.505068 -4.868567
## [31] -4.297108 -3.833194 -3.503840 -3.325433 -3.307926 -3.457655
## [37] -3.778856 -4.274127 -4.943857 -5.784457 -6.784834 -7.920367
## [43] -9.144086 -10.377572 -11.511317 -12.430793 -13.068592 -13.440575
## [49] -13.624521 -13.706896 -13.748785 -13.782319 -13.821469 -13.872769
## [55] -13.939776 -14.020358 -14.098741 -14.136324 -14.068719 -13.820057
## [61] -13.338769 -12.633433 -11.773842 -10.856577 -9.969285 -9.175272
## [67] -8.514349 -8.009317 -7.672104 -7.507927 -7.517585 -7.698307
## [73] -8.043468 -8.541337 -9.172968 -9.909643 -10.711031 -11.526432
## [79] -12.302013 -12.994333 -13.583945 -14.078851 -14.503984 -14.884222
## [85] -15.231765 -15.543263 -15.806171 -16.010199 -16.157419 -16.264982
## [91] -16.359111 -16.464519 -16.594649 -16.745743 -16.895822 -17.009646
## [97] -17.050528 -16.996280 -16.851085 -16.645161 -16.423033 -16.229091
## [103] -16.097778 -16.049668 -16.091173 -16.215626 -16.404893 -16.631990
## [109] -16.865795 -17.078322 -17.253005 -17.390290 -17.507267 -17.631541
## [115] -17.792797 -18.015496 -18.313830 -18.688042 -19.120778 -19.573264
## [121] -19.984024 -20.276769 -20.383632# Check performance
preds <- fore.aruma.wge(la$positivity, phi = e$phi, theta = e$theta, s=7, n.ahead = 7, lastn = T, limits = F)
eval_model(la$positivity,preds$f, model_name = 'AR(15) With Weekly Trend', AIC_val = e$aic) #ASE .000647
preds <- fore.aruma.wge(la$positivity, phi = e$phi, theta = e$theta, s=7, n.ahead = 90, lastn = T, limits = F)
eval_model(la$positivity,preds$f, model_name = 'AR(15) With Weekly Trend', AIC_val = e$aic) #ASE .000202
rolling <- rolling_ASE(la, e, s=7, horizon=7, model_name = 'AR(15) with Weekly Trend') #ASE .000277## [1] "test window: 240:246, train window: 210:239"## [1] "Window ASE: 0.000646772835406658"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 210:216, train window: 180:209"## [1] "Window ASE: 0.00011534695778759"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 180:186, train window: 150:179"## [1] "Window ASE: 0.0000625442246018444"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 150:156, train window: 120:149"## [1] "Window ASE: 0.000474435722736831"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 120:126, train window: 90:119"## [1] "Window ASE: 0.000099788150409411"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 90:96, train window: 60:89"## [1] "Window ASE: 0.000265503096318144"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## Warning: Removed 89 row(s) containing missing values (geom_path).Fitting additional models:
Weekly ARUMA(14,1,0)
ARUMA(14,1,0) with a weekly trend shows a slight improvement in the 7 day prediction, but a very poor performance in the 90 day trend. This is due to the 1-B term induced by the difference, that causes the trend to continue in a downward trend seen in the training data.
aic5.wge(la.d1.s7, p=0:15, q=0:2, type='bic') #BIC recommends 14,0## ---------WORKING... PLEASE WAIT...
##
##
## Five Smallest Values of bic## p q bic
## 43 14 0 -8.248609
## 46 15 0 -8.234827
## 44 14 1 -8.230948
## 45 14 2 -8.228824
## 47 15 1 -8.214616e <- est.arma.wge(la.d1.s7,p = 14, q=0)##
## Coefficients of Original polynomial:
## -0.2868 -0.0488 0.1346 -0.0076 0.0328 0.1071 -0.5340 -0.1427 -0.0600 -0.1284 -0.2537 -0.2335 -0.0086 -0.3409
##
## Factor Roots Abs Recip System Freq
## 1-1.8991B+0.9887B^2 0.9604+-0.2985i 0.9944 0.0480
## 1+0.9691B+0.9260B^2 -0.5233+-0.8978i 0.9623 0.3340
## 1+1.3945B+0.9133B^2 -0.7634+-0.7156i 0.9557 0.3801
## 1+1.8695B+0.8998B^2 -1.0388+-0.1793i 0.9486 0.4728
## 1-0.1828B+0.8139B^2 0.1123+-1.1027i 0.9022 0.2339
## 1-1.3169B+0.7598B^2 0.8665+-0.7518i 0.8717 0.1137
## 1-0.5477B+0.7325B^2 0.3738+-1.1070i 0.8558 0.1982
##
## preds <- fore.aruma.wge(la$positivity, phi = e$phi, theta=e$theta, d=1, s=7, n.ahead = 7, lastn = T, limits = F)
eval_model(la$positivity,preds$f, preds$ul, preds$ll,
'ARIMA(14,1,0) With Weekly Trend', AIC_val = e$aic) #ASE = .000541
preds <- fore.aruma.wge(la$positivity, phi = e$phi, theta=e$theta, d=1, s=7, n.ahead = 90, lastn = T, limits = F)
eval_model(la$positivity,preds$f, preds$ul, preds$ll,
'ARIMA(14,1,0) With Weekly Trend', AIC_val = e$aic) #ASE = .014017 (BAD)
rolling <- rolling_ASE(la, e, s=7, horizon=7, model_name = 'ARIMA(14,1,0) with Weekly Trend') #ASE .000679## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.00127641370441219"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 210:216, train window: 180:209"
## [1] "Window ASE: 0.000113147657616594"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 180:186, train window: 150:179"
## [1] "Window ASE: 0.000384666162028781"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 150:156, train window: 120:149"
## [1] "Window ASE: 0.00145540483664713"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 120:126, train window: 90:119"
## [1] "Window ASE: 0.000502713903767943"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 90:96, train window: 60:89"
## [1] "Window ASE: 0.000340185810696027"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## Warning: Removed 89 row(s) containing missing values (geom_path).
Signal Plus Noise
The signal plus noise model performs very well in the 7 day window, but very poorly over the 90 day window. This is a deterministic signal with a stationary mean model, and our results aren’t surprising since it doesn’t account for the weekly trend well.
########## Sig Plus Noise ##########
preds <- fore.sigplusnoise.wge(la$positivity, max.p = 15, n.ahead = 7, limits=F)##
## Coefficients of Original polynomial:
## 0.7543 0.0962 0.2365 -0.2347 -0.0441 0.1686 0.3945 -0.3197 -0.2312 -0.0513 0.0608 0.0299 0.0172 0.2035 -0.1219
##
## Factor Roots Abs Recip System Freq
## 1-1.2292B+0.9629B^2 0.6383+-0.7944i 0.9813 0.1423
## 1-0.9726B 1.0282 0.9726 0.0000
## 1+0.4746B+0.9414B^2 -0.2521+-0.9993i 0.9703 0.2893
## 1+1.6350B+0.8406B^2 -0.9725+-0.4938i 0.9168 0.4252
## 1-1.6901B+0.7766B^2 1.0882+-0.3218i 0.8812 0.0458
## 1+1.0723B+0.7471B^2 -0.7177+-0.9075i 0.8643 0.3565
## 1+0.8363B -1.1958 0.8363 0.5000
## 1-0.3085B+0.5926B^2 0.2603+-1.2727i 0.7698 0.2179
## 1-0.5721B 1.7479 0.5721 0.0000
##
## 
eval_model(la$positivity,preds$f, preds$ul, preds$ll,'SigPlusNoise', 0) #ASE = .000385
preds <- fore.sigplusnoise.wge(la$positivity, max.p = 15, n.ahead = 90, limits=F)##
## Coefficients of Original polynomial:
## 0.7576 0.0956 0.2440 -0.2442 -0.0649 0.1687 0.4054 -0.3083 -0.2337 -0.0644 0.0805 0.0514 0.0181 0.1827 -0.1313
##
## Factor Roots Abs Recip System Freq
## 1-1.2254B+0.9632B^2 0.6361+-0.7960i 0.9814 0.1427
## 1-0.9704B 1.0305 0.9704 0.0000
## 1+0.4720B+0.9372B^2 -0.2518+-1.0018i 0.9681 0.2892
## 1+1.6345B+0.8382B^2 -0.9750+-0.4924i 0.9155 0.4256
## 1-1.6494B+0.7501B^2 1.0995+-0.3527i 0.8661 0.0494
## 1+1.0698B+0.7462B^2 -0.7168+-0.9090i 0.8639 0.3563
## 1+0.8339B -1.1992 0.8339 0.5000
## 1-0.2651B+0.5828B^2 0.2275+-1.2900i 0.7634 0.2222
## 1-0.6575B 1.5209 0.6575 0.0000
##
## 
eval_model(la$positivity,preds$f, preds$ul, preds$ll,'SigPlusNoise', 0) #ASE = .003437 (BAD)
rolling <- rolling_ASE(la, e, s=7, horizon=7, model_name = 'SigPlusNoise', model_type = 'SigPlusNoise', p = 15) #ASE .000418## [1] "test window: 240:246, train window: 210:239"
##
## Coefficients of Original polynomial:
## 0.4613 -0.6072 0.2538 -0.2611 -0.0181 -0.1028 0.3613 -0.1860 -0.1092 -0.0082 -0.2988
##
## Factor Roots Abs Recip System Freq
## 1-1.2186B+0.9706B^2 0.6278+-0.7976i 0.9852 0.1439
## 1+0.4791B+0.9365B^2 -0.2558+-1.0012i 0.9677 0.2898
## 1-1.7268B+0.8266B^2 1.0445+-0.3446i 0.9092 0.0507
## 1+1.4620B+0.8123B^2 -0.8999+-0.6490i 0.9013 0.4006
## 1-0.1951B+0.6631B^2 0.1471+-1.2191i 0.8143 0.2309
## 1+0.7381B -1.3549 0.7381 0.5000
##
## 
## [1] "Window ASE: 0.000799097250179986"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 210:216, train window: 180:209"
##
## Coefficients of Original polynomial:
## 0.2158 -0.3047 -0.1716 -0.2450 -0.2585 -0.0610 0.5038 -0.3538
##
## Factor Roots Abs Recip System Freq
## 1+0.4915B+0.9561B^2 -0.2570+-0.9899i 0.9778 0.2904
## 1-1.1763B+0.9488B^2 0.6199+-0.8183i 0.9741 0.1468
## 1+1.6727B+0.8580B^2 -0.9748+-0.4640i 0.9263 0.4293
## 1-1.2037B+0.4546B^2 1.3239+-0.6686i 0.6742 0.0744
##
## 
## [1] "Window ASE: 0.000222706814850399"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 180:186, train window: 150:179"
##
## Coefficients of Original polynomial:
## 0.1201 -0.2251 -0.1223 -0.1370 -0.2003 0.2381 0.3289
##
## Factor Roots Abs Recip System Freq
## 1-1.1803B+0.9519B^2 0.6199+-0.8162i 0.9757 0.1466
## 1+0.5098B+0.8123B^2 -0.3138+-1.0642i 0.9013 0.2956
## 1-0.8070B 1.2392 0.8070 0.0000
## 1+1.3574B+0.5270B^2 -1.2878+-0.4888i 0.7260 0.4423
##
## 
## [1] "Window ASE: 0.000112207413882515"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 150:156, train window: 120:149"
##
## Coefficients of Original polynomial:
## 0.6905 -0.4220 0.1930 -0.2575 -0.0349 0.1692 0.3951 -0.3312
##
## Factor Roots Abs Recip System Freq
## 1-1.1862B+0.9627B^2 0.6161+-0.8119i 0.9812 0.1467
## 1+0.4931B+0.9174B^2 -0.2687+-1.0089i 0.9578 0.2914
## 1+1.4515B+0.6808B^2 -1.0660+-0.5766i 0.8251 0.4211
## 1-1.4488B+0.5508B^2 1.3152+-0.2928i 0.7422 0.0349
##
## 
## [1] "Window ASE: 0.000207049686158031"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 120:126, train window: 90:119"
##
## Coefficients of Original polynomial:
## 0.5644 -0.5562 0.0205 -0.1355 -0.0901 -0.0984 0.1800 -0.0990 -0.3832 0.2435 -0.3451
##
## Factor Roots Abs Recip System Freq
## 1-1.1705B+0.9263B^2 0.6318+-0.8249i 0.9624 0.1460
## 1+0.3801B+0.9120B^2 -0.2084+-1.0262i 0.9550 0.2819
## 1+1.4687B+0.8672B^2 -0.8468+-0.6604i 0.9312 0.3946
## 1-1.7165B+0.8443B^2 1.0165+-0.3888i 0.9189 0.0581
## 1+0.8842B -1.1310 0.8842 0.5000
## 1-0.4104B+0.6311B^2 0.3251+-1.2161i 0.7944 0.2084
##
## 
## [1] "Window ASE: 0.000603111475107872"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 90:96, train window: 60:89"
##
## Coefficients of Original polynomial:
## -0.1500 -0.0283 0.0601 -0.1501 -0.3244 0.0988 0.3212 -0.0211 -0.3094 -0.2922
##
## Factor Roots Abs Recip System Freq
## 1-1.1602B+0.9597B^2 0.6044+-0.8226i 0.9796 0.1491
## 1+0.4304B+0.8779B^2 -0.2451+-1.0387i 0.9370 0.2869
## 1-1.6909B+0.8047B^2 1.0506+-0.3727i 0.8970 0.0542
## 1+1.5891B+0.6861B^2 -1.1580+-0.3413i 0.8283 0.4544
## 1+0.9816B+0.6281B^2 -0.7814+-0.9907i 0.7925 0.3563
##
## 
## [1] "Window ASE: 0.00056630865776623"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## Warning: Removed 89 row(s) containing missing values (geom_path).
Neural Network
The neural network model performs fairly well, and picks up on the weekly trend, but regresses towards a mean, and is edged out by the AR(15) model from earlier.
ts_la <- ts(la$positivity[1:239], start = '1')
x = mlp(ts_la, lags = 7, hd.auto.type = 'cv', reps=10)
plot(x)
preds <- predict(x, 7)
plot(preds)
eval_model(la$positivity,preds$mean,7,model_name = 'NNFOR', AIC_val = 0) #ASE = .000766## Warning: attributes are not identical across measure variables;
## they will be dropped
ts_la <- ts(la$positivity[1:156], start = '1')
x = mlp(ts_la, lags = 7, m = 1, hd.auto.type = 'cv')
plot(x)
preds <- predict(x, 90)
plot(preds)
eval_model(la$positivity,preds$mean,90,model_name = 'NNFOR', AIC_val = 0) #ASE = .000852## Warning: attributes are not identical across measure variables;
## they will be dropped
rolling <- rolling_ASE(la, e, s=7, horizon=7, model_name = 'NNFOR', model_type = 'NNFOR') #ASE .000637## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.000982985421197367"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 210:216, train window: 180:209"
## [1] "Window ASE: 0.000125292161171138"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 180:186, train window: 150:179"
## [1] "Window ASE: 0.000157626300527327"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 150:156, train window: 120:149"
## [1] "Window ASE: 0.000171504682305925"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 120:126, train window: 90:119"
## [1] "Window ASE: 0.00158874760031172"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 90:96, train window: 60:89"
## [1] "Window ASE: 0.000544659687271713"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## Warning: Removed 89 row(s) containing missing values (geom_path).
United States
Plot the Data:
# Plot the data
plotts.sample.wge(us$positivity)
## $autplt
## [1] 1.0000000 0.9470164 0.9341278 0.9260111 0.9198625 0.9088272 0.8909389
## [8] 0.8836047 0.8738064 0.8521265 0.8325090 0.8001230 0.7724089 0.7494712
## [15] 0.7330396 0.6973509 0.6643364 0.6399202 0.6205127 0.5900024 0.5527729
## [22] 0.5218086 0.4844797 0.4575403 0.4134570 0.3899807
##
## $freq
## [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
## [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
## [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
## [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
## [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
## [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
## [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
## [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
## [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
## [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
## [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
## [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
## [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
## [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
## [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
## [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
## [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
##
## $db
## [1] 16.928178 17.041622 10.954209 3.357603 1.069112 -3.238369
## [7] -11.427364 -14.831916 -13.181967 -26.879328 -10.640049 -13.104348
## [13] -6.897898 -10.769553 -12.587656 -16.736482 -15.880304 -19.355323
## [19] -19.864320 -18.877313 -12.131779 -10.974583 -16.180375 -7.015352
## [25] -8.700809 -8.412062 -7.545175 -6.229323 -11.252124 -16.095149
## [31] -15.372565 -27.552336 -10.346212 -9.057794 -11.461366 -7.109201
## [37] -14.240281 -13.075078 -23.877820 -15.330952 -14.475847 -20.765059
## [43] -13.129351 -25.189770 -19.078147 -17.562843 -19.295449 -16.258703
## [49] -12.990975 -10.740763 -14.551855 -9.308655 -7.444421 -5.100658
## [55] -8.399431 -9.881040 -17.975943 -17.462588 -17.722631 -13.792923
## [61] -13.293199 -9.316456 -15.869391 -14.012814 -14.465016 -12.635251
## [67] -8.756403 -15.569710 -11.708858 -22.389341 -7.883153 -18.081602
## [73] -16.077320 -13.828997 -15.999068 -10.429352 -26.201640 -19.135996
## [79] -14.542536 -16.213531 -16.985225 -20.728879 -15.016682 -31.829738
## [85] -14.128733 -13.936455 -8.979219 -13.464535 -15.769438 -12.144105
## [91] -14.497446 -17.272748 -12.667286 -11.379894 -13.465174 -9.903788
## [97] -13.261966 -12.491846 -11.986603 -11.921460 -13.694377 -12.698553
## [103] -16.097000 -15.677901 -11.088504 -11.527014 -9.506795 -12.681243
## [109] -14.748674 -18.413266 -16.493124 -18.659561 -18.956181 -17.539240
## [115] -23.088025 -19.681462 -12.465460 -11.043980 -10.130884 -8.698403
## [121] -17.478083 -17.827845 -23.241743
##
## $dbz
## [1] 12.8975742 12.6104343 12.1288216 11.4480751 10.5614946 9.4601976
## [7] 8.1330155 6.5665777 4.7459620 2.6568532 0.2916097 -2.3347667
## [13] -5.1453437 -7.9215132 -10.2286515 -11.5956090 -11.9705143 -11.6983430
## [19] -11.1158574 -10.4346224 -9.7877324 -9.2573489 -8.8845067 -8.6785829
## [25] -8.6271903 -8.7045894 -8.8785406 -9.1161191 -9.3884024 -9.6733156
## [31] -9.9562335 -10.2290872 -10.4895332 -10.7412048 -10.9945308 -11.2664674
## [37] -11.5777577 -11.9475163 -12.3857167 -12.8839971 -13.4053244 -13.8761561
## [43] -14.1925522 -14.2554748 -14.0264019 -13.5536343 -12.9396293 -12.2881747
## [49] -11.6760017 -11.1508180 -10.7394607 -10.4558243 -10.3057268 -10.2889439
## [55] -10.3991649 -10.6224555 -10.9347537 -11.2993911 -11.6668145 -11.9799663
## [61] -12.1877268 -12.2626791 -12.2120173 -12.0727873 -11.8955891 -11.7286297
## [67] -11.6092436 -11.5622931 -11.6019261 -11.7339670 -11.9576927 -12.2666738
## [73] -12.6487555 -13.0854284 -13.5510287 -14.0125543 -14.4313384 -14.7679001
## [79] -14.9901220 -15.0821781 -15.0493866 -14.9156322 -14.7148182 -14.4811756
## [85] -14.2423645 -14.0164221 -13.8117305 -13.6288608 -13.4634734 -13.3096700
## [91] -13.1631127 -13.0231056 -12.8930232 -12.7790349 -12.6877269 -12.6235629
## [97] -12.5870636 -12.5742955 -12.5779435 -12.5898586 -12.6044544 -12.6217762
## [103] -12.6489228 -12.6990700 -12.7883489 -12.9315227 -13.1373601 -13.4041478
## [109] -13.7156159 -14.0383220 -14.3232947 -14.5159238 -14.5746988 -14.4899695
## [115] -14.2884211 -14.0193071 -13.7340281 -13.4719928 -13.2561370 -13.0948841
## [121] -12.9868227 -12.9259627 -12.9064607# Look deeper at spectral density
parzen.wge(us$positivity, trun=100)## $freq
## [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
## [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
## [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
## [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
## [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
## [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
## [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
## [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
## [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
## [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
## [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
## [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
## [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
## [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
## [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
## [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
## [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
##
## $pzgram
## [1] 15.382933 14.588356 12.831036 9.869648 5.787060 1.193729
## [7] -3.198010 -7.039404 -9.463870 -10.954884 -11.649927 -11.146329
## [13] -10.497311 -10.387322 -10.996468 -12.241666 -14.051773 -16.357297
## [19] -17.782568 -16.624195 -14.478599 -12.431116 -10.590963 -9.053923
## [25] -7.897040 -7.190547 -7.099666 -7.820660 -9.477927 -12.016628
## [31] -14.502154 -14.318422 -11.893364 -9.673811 -8.566444 -8.883526
## [37] -10.680694 -13.468736 -15.662532 -16.381447 -16.513312 -16.415800
## [43] -16.292762 -16.362416 -16.490077 -16.180717 -15.264532 -14.173466
## [49] -13.172032 -12.002961 -10.387418 -8.753405 -7.762275 -7.720243
## [55] -8.641456 -10.426594 -12.887273 -15.145752 -15.232845 -13.814286
## [61] -12.874594 -12.869708 -13.340726 -13.353437 -12.555602 -11.409945
## [67] -10.418409 -9.928099 -9.987950 -10.427637 -11.256579 -12.668644
## [73] -14.277216 -14.778950 -14.334108 -14.338547 -14.997710 -15.512866
## [79] -15.646402 -16.183486 -17.184966 -17.379811 -16.333944 -15.008996
## [85] -13.685617 -12.594381 -12.257052 -12.870686 -13.936525 -14.509900
## [91] -14.398784 -13.976794 -13.305756 -12.477274 -11.806977 -11.561958
## [97] -11.724854 -12.038907 -12.361702 -12.843276 -13.613715 -14.408198
## [103] -14.402867 -13.168262 -11.604517 -10.672442 -10.804908 -12.107197
## [109] -14.317964 -16.503509 -17.646159 -18.062309 -18.410008 -18.446313
## [115] -17.213576 -15.075943 -13.043248 -11.492666 -10.719546 -11.075606
## [121] -12.778110 -15.616148 -17.533541
There is surpsingly little evidence of serial correlation in the US data compared to Louisiana. There are very weak peaks in the spectral density plots and a very slowly damping behavior in the autocorrelation function plots. Increasing the truncation point of the Parzen function doesn’t show significant points to indicate a significant serial correlation.
Evaluate Stationarity:
There are few indications that suggest the data is not stationary, it appears that this could be a stationary time series with some wandering behavior.
Model Construction:
Differencing the data once shows no significant cyclical pattern in the ACF.
Identify correlation structures:
Here we are going to look at different filters to try to model out some of the correlation in the data. We’ll also look at a factor table to see if we can identify a weekly or monthly trend.
# Difference the data once
us.d1 <- artrans.wge(us$positivity,1)
# Difference it again, no improvement
artrans.wge(us.d1,1)
## [1] -0.01425418201 0.01978404489 0.00344009601 -0.06802340221 0.07352438226
## [6] 0.00692990610 -0.00714591876 -0.02084688094 -0.00243451668 0.04964910774
## [11] -0.03317545244 -0.05925122282 0.04131133599 0.02405539404 0.00601164353
## [16] -0.01817377270 0.03083149410 -0.08610403374 0.08409039403 -0.02307959862
## [21] -0.01642926846 0.01098569962 -0.09811599954 0.15692527231 -0.08792327276
## [26] 0.03274987079 -0.01386236757 0.00918935394 0.00950420531 -0.01219184870
## [31] -0.02247301575 0.00493196641 -0.01035396435 0.06647902814 -0.06646866380
## [36] 0.04033066349 -0.03339489118 0.02054824224 -0.00647801193 -0.01504550299
## [41] -0.05322120943 0.13931525485 -0.07895743252 -0.00049040986 0.01314512699
## [46] -0.01923457995 0.02586767756 -0.02425491425 0.02450886876 -0.00465228139
## [51] -0.01169959290 0.01236632704 -0.02144306358 0.00542332225 0.02043214764
## [56] -0.01299864507 0.00739442008 -0.00979988508 0.00423135860 -0.02960535336
## [61] 0.05470381881 -0.02719446153 0.00982324167 -0.01199966873 0.00805208606
## [66] -0.00623128961 0.00376016574 -0.00321244605 0.00046827232 0.00967439924
## [71] -0.00039847540 -0.01103591006 0.00557147233 -0.00879606612 0.01939398534
## [76] -0.01744094841 0.00899981221 -0.00453842682 0.00960374829 -0.00902539121
## [81] 0.00071623524 -0.00215174346 0.00057748709 0.00145547567 -0.00053317682
## [86] 0.00332735306 0.00265811119 -0.00623099656 -0.00110895510 0.01109982666
## [91] -0.00786465370 -0.00713998322 0.01679738401 -0.01340796604 0.00140554934
## [96] 0.00886218895 -0.00864097890 0.00711356390 -0.00487214075 0.00488677486
## [101] -0.00584010118 0.00130517473 -0.00229574827 0.02632837457 -0.04012080831
## [106] 0.01463289507 0.00921421538 -0.00786276163 -0.00956919410 0.02247511852
## [111] -0.00867981153 -0.00477364349 -0.00200677029 0.01103878341 -0.02503295136
## [116] 0.01482554512 0.01315600697 0.00117340963 -0.02472644128 0.01051180247
## [121] 0.00847784399 -0.01529062714 0.01004135529 -0.00054810898 -0.00068255160
## [126] -0.00003667196 0.00499935746 -0.01283051447 0.00653668209 0.00001297137
## [131] 0.00435733193 -0.00366747653 0.00009589148 -0.00386304015 0.00230426213
## [136] -0.00426560014 0.00520113038 0.00308756001 -0.00116121005 0.00533521705
## [141] -0.01614762899 0.01470606243 -0.01621035783 0.01068486672 0.00547535911
## [146] -0.00548022191 -0.00180068959 0.00571874086 -0.00010420098 -0.01070456362
## [151] -0.00465714413 0.02870969236 0.01514774804 -0.08250772757 0.05829342260
## [156] -0.00198719651 -0.02080531661 0.01639823857 0.00370790673 -0.00407842935
## [161] -0.00434054222 0.00265285183 -0.00085254569 0.00250766698 -0.00625128498
## [166] 0.00942351027 0.00338309003 -0.01553426756 0.00769235953 -0.00118990894
## [171] -0.00037389962 -0.00315165280 0.01446704665 -0.02152086746 0.02277328413
## [176] -0.01200183966 -0.00177860185 -0.00346617760 0.00911724230 -0.00659330821
## [181] 0.01005544004 -0.00283078439 -0.00248659228 -0.00708970815 0.00880266567
## [186] 0.02081367549 -0.04811764792 0.03040237498 -0.00930103949 0.00235031086
## [191] -0.00360694153 0.00316510150 0.01613333560 -0.00720811766 -0.02667476139
## [196] 0.02312455662 0.00430214086 -0.01507066638 0.00371924998 -0.00424231109
## [201] 0.01639379996 0.00058917864 -0.01267489139 -0.00217868997 0.01234496297
## [206] -0.01613330064 0.01052661225 0.00344839781 0.00460398622 -0.01318871581
## [211] 0.00352505023 -0.00084230457 -0.00159297749 0.00027587920 0.00862199954
## [216] 0.00155555540 -0.00407252270 -0.00513442594 -0.00960502040 0.01307353256
## [221] 0.00264711510 0.00880620228 -0.01204944398 -0.00666575615 0.00578262594
## [226] 0.00311978586 -0.01173173304 0.00533528809 0.01865195916 -0.01111606447
## [231] -0.00982475902 0.01001397542 -0.01088472895 0.00645926430 -0.00394642266
## [236] 0.01856996819 -0.00723121433 -0.01788513476 0.01628207270 0.00493381307
## [241] -0.01441135849 -0.00594183904 0.03534821567 -0.02865985992# Try various seasonal models
us.s7 <- artrans.wge(us$positivity, c(rep(0,6),1))
us.s12 <- artrans.wge(us$positivity, c(rep(0,11),1))
us.d1.s7 <- artrans.wge(us.d1, c(rep(0,6),1))
us.d1.s12 <- artrans.wge(us.d1, c(rep(0,11),1))
# Overfit table
est.ar.wge(us$positivity, p=15)##
## Coefficients of Original polynomial:
## 0.3866 0.1849 0.1179 0.1508 0.1336 -0.0545 0.1364 0.2048 0.0220 0.0994 -0.1207 -0.1809 -0.0703 0.1205 -0.1875
##
## Factor Roots Abs Recip System Freq
## 1-1.9642B+0.9671B^2 1.0155+-0.0524i 0.9834 0.0082
## 1-0.3434B+0.8525B^2 0.2014+-1.0642i 0.9233 0.2202
## 1+0.9180B -1.0893 0.9180 0.5000
## 1+0.4309B+0.8427B^2 -0.2557+-1.0589i 0.9180 0.2877
## 1-1.3843B+0.8248B^2 0.8392+-0.7129i 0.9082 0.1121
## 1+1.2323B+0.8146B^2 -0.7564+-0.8096i 0.9025 0.3696
## 1+1.6110B+0.8080B^2 -0.9970+-0.4937i 0.8989 0.4268
## 1-0.8868B+0.5415B^2 0.8189+-1.0845i 0.7358 0.1471
##
##
##
## Coefficients of Original polynomial:
## 0.3866 0.1849 0.1179 0.1508 0.1336 -0.0545 0.1364 0.2048 0.0220 0.0994 -0.1207 -0.1809 -0.0703 0.1205 -0.1875
##
## Factor Roots Abs Recip System Freq
## 1-1.9642B+0.9671B^2 1.0155+-0.0524i 0.9834 0.0082
## 1-0.3434B+0.8525B^2 0.2014+-1.0642i 0.9233 0.2202
## 1+0.9180B -1.0893 0.9180 0.5000
## 1+0.4309B+0.8427B^2 -0.2557+-1.0589i 0.9180 0.2877
## 1-1.3843B+0.8248B^2 0.8392+-0.7129i 0.9082 0.1121
## 1+1.2323B+0.8146B^2 -0.7564+-0.8096i 0.9025 0.3696
## 1+1.6110B+0.8080B^2 -0.9970+-0.4937i 0.8989 0.4268
## 1-0.8868B+0.5415B^2 0.8189+-1.0845i 0.7358 0.1471
##
## ## $phi
## [1] 0.38657693 0.18487981 0.11787070 0.15076990 0.13357768 -0.05449747
## [7] 0.13638095 0.20481227 0.02200580 0.09936180 -0.12065274 -0.18089166
## [13] -0.07032510 0.12048686 -0.18746678
##
## $res
## [1] 0.01658565790 0.01375283230 -0.00090975462 0.00845540202 0.01787093204
## [6] -0.03984666673 -0.00642488247 0.01431184179 0.01779142314 0.01080242259
## [11] 0.00663355927 0.04102041625 0.03926094547 -0.01183329929 -0.01332829485
## [16] 0.01014091369 0.01845543645 0.00884849996 0.04057279959 -0.02456491855
## [21] 0.01134078334 0.01787148894 0.00870671436 0.02431789715 -0.07283238027
## [26] -0.00517819506 -0.03437950208 0.00118715503 -0.01032289425 0.00644319619
## [31] 0.00652596746 0.01987781181 0.00488724220 0.00373016664 -0.01721744312
## [36] 0.02418732184 -0.01049618328 0.01171728082 0.00293030831 -0.01518432094
## [41] -0.00155289996 -0.01039100053 -0.06795034138 0.02871827380 0.00452748430
## [46] -0.01239578134 0.00854411475 -0.00827701124 0.00216160779 -0.00435294323
## [51] 0.02690416243 0.01475487501 0.01892083537 0.00496840193 -0.01299866922
## [56] -0.00935331316 0.01502989285 -0.00836566651 0.00641615521 0.00149736569
## [61] -0.00310893443 -0.03220818092 0.00693531442 0.00508118048 0.00916087672
## [66] 0.00403499757 0.00518738405 -0.00595027914 0.00420244633 0.00369506342
## [71] -0.00360685454 0.00856749870 0.00556391548 -0.00256806545 -0.00069150154
## [76] -0.00312190783 0.00158672696 -0.00127510904 0.00152979208 -0.00120194631
## [81] 0.00433291352 -0.00098428801 -0.00156131332 -0.00254140231 -0.00486746293
## [86] -0.00545535807 -0.00710667225 -0.00227452103 0.00061060139 -0.00133946919
## [91] -0.00769660571 0.00591309232 0.00306698063 -0.00369564838 0.00718262066
## [96] -0.00001671539 -0.00461294328 0.00401196919 -0.00159613819 0.00109499152
## [101] 0.00196209263 0.00371501321 -0.00167261727 0.00186236438 -0.00299427170
## [106] 0.01941514673 -0.00382930045 -0.00384991841 0.00190630333 -0.00071041088
## [111] -0.01224783761 0.00632026187 0.00502788605 -0.00259796971 -0.00119884956
## [116] 0.00344108999 -0.01517812106 -0.00822784957 0.00497434318 0.00772444255
## [121] -0.00229308135 -0.00430162836 -0.00085720848 -0.00944426026 0.00157793679
## [126] 0.00018534116 -0.00056310649 -0.00126815864 0.00390301921 -0.00977286373
## [131] -0.00092462071 -0.00053444002 -0.00068972220 -0.00069105253 0.00380122262
## [136] -0.00448148251 -0.00643582696 -0.00728575132 -0.00634107135 -0.00040158021
## [141] 0.00118932466 0.00684192504 -0.00612621632 0.00392189653 -0.00818374001
## [146] -0.00325862021 0.00308286413 0.00267336077 -0.00223351023 0.00232775916
## [151] 0.00372790022 -0.00494997796 -0.01113158303 0.00976643002 0.03693895275
## [156] -0.02817729817 -0.00566048758 -0.00240558087 -0.01928741546 -0.00913914196
## [161] 0.00412288863 -0.00488766390 -0.00643902903 -0.00091385188 -0.00827121966
## [166] 0.00636653750 0.00466488156 0.00003018588 0.00509632278 0.00592229541
## [171] -0.00489662044 -0.00466327616 -0.00237286795 -0.00815518027 0.00495705805
## [176] -0.01135691942 0.00368596891 0.00028501350 -0.00433408475 -0.00817656855
## [181] 0.00080322115 -0.00704205020 0.00114529007 0.00694825730 0.00090665491
## [186] -0.00526555278 -0.00078038882 0.02232359024 -0.01211795724 0.00344900910
## [191] -0.00690745033 -0.00855145879 -0.01250710823 -0.00693007971 0.00369252631
## [196] 0.00846396477 -0.01386141770 -0.00732336885 0.00637722556 -0.00252184851
## [201] -0.00431535013 -0.01269272632 0.00126723235 0.00387963806 0.00032530003
## [206] -0.00700880871 0.00587812304 -0.00701177611 -0.00660973973 0.00352045916
## [211] 0.01335070103 -0.00109421533 -0.00276350971 -0.00446999423 -0.00667799040
## [216] -0.00394059592 0.00040448468 0.00528969534 0.00592415760 0.00070268984
## [221] -0.01447748759 -0.00501264793 0.00045067434 0.01262251779 0.00709674750
## [226] -0.00096957046 -0.00472347849 -0.00187216773 -0.00895473279 -0.00658333715
## [231] 0.01023025936 0.00504505014 -0.00489649435 0.00159401036 -0.00531932415
## [236] -0.00433809527 -0.00418639105 0.00973786582 0.01178895197 -0.00241604198
## [241] 0.00069451250 0.00853576955 0.00406154919 -0.00520709593 0.02140959751
## [246] 0.00587643257
##
## $avar
## [1] 0.0001537157
##
## $aic
## [1] -8.650324
##
## $aicc
## [1] -7.631283
##
## $bic
## [1] -8.422335# Generate a factor table for a (1-B^7)
tswge::factor.wge(c(rep(0,6),1))##
## Coefficients of Original polynomial:
## 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000
##
## Factor Roots Abs Recip System Freq
## 1-1.0000B 1.0000 1.0000 0.0000
## 1+0.4450B+1.0000B^2 -0.2225+-0.9749i 1.0000 0.2857
## 1-1.2470B+1.0000B^2 0.6235+-0.7818i 1.0000 0.1429
## 1+1.8019B+1.0000B^2 -0.9010+-0.4339i 1.0000 0.4286
##
## Estimate parameters and check residuals:
The ACF and PACF plots look like white noise, as expected. The Ljung-Box test fails to reject the null hyphothesis that the residuals are white noise. The generated spectral densities look like they have a fairly good fit, the final peak is overpronounced in our generations compared to the actual data. The ACF comparisons are all over the board, but with peaks in the same general spots (lags at 7, 14, 21). The generated realizations look quite a bit different, but the spectral densities and ACF plots look fairly similar, indicating this could be a useful model.
# Estimate the model params
aic5.wge(us$positivity, p=0:20, q=0:2, type='bic') # aic = 19,1 > 18,0 : bic = 1,1 > 9,1 > 2,1## ---------WORKING... PLEASE WAIT...
##
##
## Error in aic calculation at 20 1
## Five Smallest Values of bic## p q bic
## 5 1 1 -8.487898
## 29 9 1 -8.468857
## 8 2 1 -8.465559
## 6 1 2 -8.465077
## 35 11 1 -8.462634e <- est.arma.wge(us$positivity, p=9, q=1)##
## Coefficients of Original polynomial:
## 1.2711 -0.1631 -0.0699 0.0662 0.0028 -0.1636 0.1612 0.1144 -0.2295
##
## Factor Roots Abs Recip System Freq
## 1-1.9861B+0.9886B^2 1.0045+-0.0501i 0.9943 0.0079
## 1-0.0326B+0.7214B^2 0.0226+-1.1772i 0.8493 0.2469
## 1-1.1925B+0.6943B^2 0.8588+-0.8384i 0.8332 0.1231
## 1+1.2068B+0.6319B^2 -0.9548+-0.8190i 0.7949 0.3872
## 1+0.7333B -1.3636 0.7333 0.5000
##
## # Ljung Box Test still shows white noise residuals
ljung.wge(e$res)## Obs -0.03356559 -0.04604162 -0.004168052 0.01547628 -0.01688291 -0.03271148 -0.0009237203 0.003706762 0.06424109 0.1178936 -0.03902071 -0.1007295 -0.06082771 0.1679569 -0.04072279 -0.1047464 -0.02937219 0.1135898 0.1318257 0.02116977 0.05601096 -0.05664106 0.07570874 -0.1847351## $test
## [1] "Ljung-Box test"
##
## $K
## [1] 24
##
## $chi.square
## [1] 41.79917
##
## $df
## [1] 24
##
## $pval
## [1] 0.0135891ljung.wge(e$res, K = 48)## Obs -0.03356559 -0.04604162 -0.004168052 0.01547628 -0.01688291 -0.03271148 -0.0009237203 0.003706762 0.06424109 0.1178936 -0.03902071 -0.1007295 -0.06082771 0.1679569 -0.04072279 -0.1047464 -0.02937219 0.1135898 0.1318257 0.02116977 0.05601096 -0.05664106 0.07570874 -0.1847351 -0.04456868 -0.1097331 -0.008943643 0.1029134 0.03469855 -0.1181843 -0.04739403 0.0068443 0.01577266 -0.04474552 0.08704028 0.02497461 0.1092873 -0.08576074 0.01382154 -0.01940963 0.0472867 0.06192417 -0.01025361 -0.02817321 -0.0161036 0.04494596 -0.02163542 0.009431244## $test
## [1] "Ljung-Box test"
##
## $K
## [1] 48
##
## $chi.square
## [1] 65.38729
##
## $df
## [1] 48
##
## $pval
## [1] 0.04816801acf(e$res)
pacf(e$res)
dev.off()## null device
## 1#Compare Spectral Densities
sims = 5
SpecDen = parzen.wge(us$positivity, plot = "FALSE")
plot(SpecDen$freq, SpecDen$pzgram, type = "l", lwd = 6)
for( i in 1: sims)
{
SpecDen2 = parzen.wge(gen.aruma.wge(246,s = 0, phi = e$phi, plot ="FALSE"), plot = "FALSE")
lines(SpecDen2$freq,SpecDen2$pzgram, lwd = 2, col = "red")
}
#Compare ACFs
sims = 5
ACF = acf(us$positivity, plot = "FALSE")
plot(ACF$lag ,ACF$acf , type = "l", lwd = 6)
for( i in 1: sims)
{
ACF2 = acf(gen.aruma.wge(246, s = 0, phi = e$phi, plot = "FALSE"), plot = "FALSE")
lines(ACF2$lag ,ACF2$acf, lwd = 2, col = "red")
}
#Compare Generated Realizations
eGen = gen.aruma.wge(246, s = 0, phi = e$phi, vara = e$avar)
plotts.sample.wge(eGen)## $autplt
## [1] 1.0000000 0.9935339 0.9838052 0.9713461 0.9562024 0.9380209 0.9171481
## [8] 0.8939896 0.8687276 0.8414615 0.8122268 0.7814705 0.7491547 0.7152164
## [15] 0.6803175 0.6438795 0.6061015 0.5667866 0.5267571 0.4858514 0.4440230
## [22] 0.4014926 0.3584973 0.3151003 0.2707877 0.2260548
##
## $freq
## [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
## [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
## [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
## [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
## [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
## [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
## [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
## [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
## [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
## [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
## [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
## [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
## [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
## [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
## [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
## [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
## [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
##
## $db
## [1] 4.785980 20.460683 6.878691 -9.792577 -2.442456 -6.769525
## [7] -2.114508 -5.513326 -4.666739 -8.681357 -6.524804 -7.989461
## [13] -9.545319 -10.547635 -9.767938 -13.914741 -13.222626 -11.588570
## [19] -14.143523 -12.485914 -14.992938 -14.761362 -19.147332 -12.042511
## [25] -18.455932 -16.026129 -15.447145 -20.189811 -18.951746 -35.972499
## [31] -15.036756 -21.535024 -14.861839 -22.526411 -18.001167 -19.586483
## [37] -18.745639 -18.769383 -22.678378 -21.667724 -19.555406 -18.744277
## [43] -19.948553 -19.908577 -21.884591 -23.542054 -21.372526 -20.604435
## [49] -18.870470 -21.058358 -18.361203 -18.387750 -22.869790 -22.423839
## [55] -21.306296 -22.262206 -24.710437 -24.561472 -19.060970 -26.298557
## [61] -22.948965 -22.721761 -21.603496 -20.555001 -22.313456 -21.712554
## [67] -26.106076 -25.014777 -23.047509 -23.253366 -24.970398 -28.859453
## [73] -24.665586 -23.240980 -31.541241 -23.468987 -25.489892 -24.517900
## [79] -25.702613 -28.213259 -25.260366 -28.133953 -25.136608 -23.862811
## [85] -22.097075 -27.666242 -24.031868 -23.117571 -27.022319 -27.565200
## [91] -27.458709 -23.068928 -27.052203 -27.098103 -31.964080 -26.061005
## [97] -28.751894 -26.069454 -26.578042 -27.960504 -29.763770 -26.343833
## [103] -26.969098 -23.194316 -26.482595 -25.923156 -25.220976 -31.365838
## [109] -26.645385 -23.904860 -26.940378 -28.862589 -28.254769 -27.305370
## [115] -30.382414 -24.531366 -28.769833 -23.600321 -28.620172 -26.270184
## [121] -26.089719 -23.583943 -26.787258
##
## $dbz
## [1] 12.925445 12.659061 12.212530 11.581967 10.761845 9.744968
## [7] 8.522543 7.084561 5.420919 3.524292 1.397019 -0.933425
## [13] -3.383545 -5.769855 -7.805983 -9.251522 -10.119925 -10.619053
## [19] -10.932032 -11.147277 -11.308923 -11.465576 -11.679084 -12.009401
## [25] -12.498549 -13.161416 -13.979962 -14.897436 -15.816880 -16.618552
## [31] -17.207066 -17.561524 -17.736164 -17.811584 -17.850539 -17.886983
## [37] -17.936897 -18.012207 -18.127776 -18.299997 -18.540416 -18.848757
## [43] -19.208424 -19.586555 -19.940482 -20.230689 -20.435419 -20.557951
## [49] -20.621475 -20.656312 -20.688648 -20.735738 -20.806683 -20.905491
## [55] -21.033596 -21.190547 -21.373111 -21.573992 -21.781707 -21.982598
## [61] -22.164812 -22.322612 -22.458625 -22.582582 -22.707176 -22.843271
## [67] -22.996590 -23.166879 -23.349321 -23.537228 -23.724637 -23.907666
## [73] -24.084143 -24.252047 -24.407932 -24.546465 -24.661534 -24.748438
## [79] -24.806027 -24.837644 -24.850502 -24.853978 -24.857791 -24.870730
## [85] -24.900084 -24.951444 -25.028498 -25.132602 -25.262198 -25.412357
## [91] -25.574816 -25.738750 -25.892239 -26.024081 -26.125357 -26.190261
## [97] -26.216163 -26.203322 -26.154771 -26.076539 -25.977872 -25.870918
## [103] -25.769584 -25.687730 -25.637172 -25.625990 -25.657378 -25.729124
## [109] -25.833696 -25.958940 -26.089447 -26.208624 -26.301300 -26.356388
## [115] -26.368892 -26.340592 -26.279221 -26.196519 -26.105888 -26.020309
## [121] -25.950849 -25.905790 -25.890203plotts.sample.wge(us$positivity)## $autplt
## [1] 1.0000000 0.9470164 0.9341278 0.9260111 0.9198625 0.9088272 0.8909389
## [8] 0.8836047 0.8738064 0.8521265 0.8325090 0.8001230 0.7724089 0.7494712
## [15] 0.7330396 0.6973509 0.6643364 0.6399202 0.6205127 0.5900024 0.5527729
## [22] 0.5218086 0.4844797 0.4575403 0.4134570 0.3899807
##
## $freq
## [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
## [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
## [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
## [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
## [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
## [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
## [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
## [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
## [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
## [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
## [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
## [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
## [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
## [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
## [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
## [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
## [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
##
## $db
## [1] 16.928178 17.041622 10.954209 3.357603 1.069112 -3.238369
## [7] -11.427364 -14.831916 -13.181967 -26.879328 -10.640049 -13.104348
## [13] -6.897898 -10.769553 -12.587656 -16.736482 -15.880304 -19.355323
## [19] -19.864320 -18.877313 -12.131779 -10.974583 -16.180375 -7.015352
## [25] -8.700809 -8.412062 -7.545175 -6.229323 -11.252124 -16.095149
## [31] -15.372565 -27.552336 -10.346212 -9.057794 -11.461366 -7.109201
## [37] -14.240281 -13.075078 -23.877820 -15.330952 -14.475847 -20.765059
## [43] -13.129351 -25.189770 -19.078147 -17.562843 -19.295449 -16.258703
## [49] -12.990975 -10.740763 -14.551855 -9.308655 -7.444421 -5.100658
## [55] -8.399431 -9.881040 -17.975943 -17.462588 -17.722631 -13.792923
## [61] -13.293199 -9.316456 -15.869391 -14.012814 -14.465016 -12.635251
## [67] -8.756403 -15.569710 -11.708858 -22.389341 -7.883153 -18.081602
## [73] -16.077320 -13.828997 -15.999068 -10.429352 -26.201640 -19.135996
## [79] -14.542536 -16.213531 -16.985225 -20.728879 -15.016682 -31.829738
## [85] -14.128733 -13.936455 -8.979219 -13.464535 -15.769438 -12.144105
## [91] -14.497446 -17.272748 -12.667286 -11.379894 -13.465174 -9.903788
## [97] -13.261966 -12.491846 -11.986603 -11.921460 -13.694377 -12.698553
## [103] -16.097000 -15.677901 -11.088504 -11.527014 -9.506795 -12.681243
## [109] -14.748674 -18.413266 -16.493124 -18.659561 -18.956181 -17.539240
## [115] -23.088025 -19.681462 -12.465460 -11.043980 -10.130884 -8.698403
## [121] -17.478083 -17.827845 -23.241743
##
## $dbz
## [1] 12.8975742 12.6104343 12.1288216 11.4480751 10.5614946 9.4601976
## [7] 8.1330155 6.5665777 4.7459620 2.6568532 0.2916097 -2.3347667
## [13] -5.1453437 -7.9215132 -10.2286515 -11.5956090 -11.9705143 -11.6983430
## [19] -11.1158574 -10.4346224 -9.7877324 -9.2573489 -8.8845067 -8.6785829
## [25] -8.6271903 -8.7045894 -8.8785406 -9.1161191 -9.3884024 -9.6733156
## [31] -9.9562335 -10.2290872 -10.4895332 -10.7412048 -10.9945308 -11.2664674
## [37] -11.5777577 -11.9475163 -12.3857167 -12.8839971 -13.4053244 -13.8761561
## [43] -14.1925522 -14.2554748 -14.0264019 -13.5536343 -12.9396293 -12.2881747
## [49] -11.6760017 -11.1508180 -10.7394607 -10.4558243 -10.3057268 -10.2889439
## [55] -10.3991649 -10.6224555 -10.9347537 -11.2993911 -11.6668145 -11.9799663
## [61] -12.1877268 -12.2626791 -12.2120173 -12.0727873 -11.8955891 -11.7286297
## [67] -11.6092436 -11.5622931 -11.6019261 -11.7339670 -11.9576927 -12.2666738
## [73] -12.6487555 -13.0854284 -13.5510287 -14.0125543 -14.4313384 -14.7679001
## [79] -14.9901220 -15.0821781 -15.0493866 -14.9156322 -14.7148182 -14.4811756
## [85] -14.2423645 -14.0164221 -13.8117305 -13.6288608 -13.4634734 -13.3096700
## [91] -13.1631127 -13.0231056 -12.8930232 -12.7790349 -12.6877269 -12.6235629
## [97] -12.5870636 -12.5742955 -12.5779435 -12.5898586 -12.6044544 -12.6217762
## [103] -12.6489228 -12.6990700 -12.7883489 -12.9315227 -13.1373601 -13.4041478
## [109] -13.7156159 -14.0383220 -14.3232947 -14.5159238 -14.5746988 -14.4899695
## [115] -14.2884211 -14.0193071 -13.7340281 -13.4719928 -13.2561370 -13.0948841
## [121] -12.9868227 -12.9259627 -12.9064607# Check performance
preds <- fore.aruma.wge(us$positivity, phi = e$phi, theta = e$theta, n.ahead = 7, lastn = T, limits = F)
eval_model(us$positivity,preds$f, model_name = 'ARMA(9,1)',AIC_val = e$aic) #ASE .000121
preds <- fore.aruma.wge(us$positivity, phi = e$phi, theta = e$theta, n.ahead = 90, lastn = T, limits = F)
eval_model(us$positivity,preds$f, model_name = 'ARMA(9,1)', AIC_val = e$aic) #ASE .00115
rolling <- rolling_ASE(us, e, s=0, horizon=7, model_name = 'AR(9,1)' ) #ASE .000096## [1] "test window: 240:246, train window: 210:239"## [1] "Window ASE: 0.00019584690928318"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 210:216, train window: 180:209"## [1] "Window ASE: 0.0000286092864906074"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 180:186, train window: 150:179"## [1] "Window ASE: 0.0000184254204829249"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 150:156, train window: 120:149"## [1] "Window ASE: 0.000264139426854921"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 120:126, train window: 90:119"## [1] "Window ASE: 0.0000280660276168899"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 90:96, train window: 60:89"## [1] "Window ASE: 0.0000435508295999369"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## Warning: Removed 89 row(s) containing missing values (geom_path).Fitting additional models:
ARMA(18,1) with a weekly trend
The ARMA(18,1) with a weekly trend shows a slight drop in performance in the 7 day prediction, but much better performance in the 90 day trend. It is the most balanced model, having the best performance in 90 day predictions and competitive performance in the rolling window ASE.
aic5.wge(us.s7, p=0:20, q=0:5, type='bic') #aic = 20,1 > 20,2 > 18,1 : bic = 20,1 > 18,1## ---------WORKING... PLEASE WAIT...
##
##
## Error in aic calculation at 13 1
## Error in aic calculation at 15 1
## Error in aic calculation at 15 4
## Error in aic calculation at 15 5
## Error in aic calculation at 16 3
## Five Smallest Values of bic## p q bic
## 122 20 1 -8.432983
## 110 18 1 -8.421151
## 116 19 1 -8.403711
## 123 20 2 -8.397647
## 53 8 4 -8.354801e <- est.arma.wge(us.s7, p=18, q=1)##
## Coefficients of Original polynomial:
## 1.2038 -0.1377 0.0970 -0.1074 -0.0016 -0.0903 -0.6689 1.0683 -0.2391 0.2188 -0.4459 -0.0299 0.0073 -0.1330 0.4938 -0.2741 0.3459 -0.3470
##
## Factor Roots Abs Recip System Freq
## 1-1.9871B+0.9916B^2 1.0019+-0.0674i 0.9958 0.0107
## 1-1.5531B+0.9722B^2 0.7988+-0.6249i 0.9860 0.1057
## 1+1.9252B+0.9394B^2 -1.0247+-0.1205i 0.9692 0.4814
## 1+1.1191B+0.9092B^2 -0.6154+-0.8492i 0.9535 0.3498
## 1+1.4737B+0.9061B^2 -0.8132+-0.6651i 0.9519 0.3909
## 1-1.7816B+0.9039B^2 0.9855+-0.3675i 0.9507 0.0568
## 1-0.3277B+0.8822B^2 0.1857+-1.0484i 0.9392 0.2221
## 1-0.5258B+0.8490B^2 0.3096+-1.0402i 0.9214 0.2040
## 1+0.4535B+0.6870B^2 -0.3301+-1.1604i 0.8289 0.2941
##
## # Ljung Box Test still shows white noise residuals
ljung.wge(artrans.wge(us.s7, phi.tr = e$phi))## Obs -0.5065971 -0.003026183 0.002316859 0.08036559 -0.1209084 0.07023771 0.04574104 -0.03331472 0.03382918 -0.103711 0.001563584 0.07111938 0.1553696 -0.2417995 -0.04760502 0.2130846 -0.1512536 0.01712802 0.1140875 0.00473265 -0.1590841 0.1094058 -0.09090864 0.1008303## $test
## [1] "Ljung-Box test"
##
## $K
## [1] 24
##
## $chi.square
## [1] 121.815
##
## $df
## [1] 24
##
## $pval
## [1] 0.000000000000004551914ljung.wge(artrans.wge(us.s7, phi.tr = e$phi), K = 48)
## Obs -0.5065971 -0.003026183 0.002316859 0.08036559 -0.1209084 0.07023771 0.04574104 -0.03331472 0.03382918 -0.103711 0.001563584 0.07111938 0.1553696 -0.2417995 -0.04760502 0.2130846 -0.1512536 0.01712802 0.1140875 0.00473265 -0.1590841 0.1094058 -0.09090864 0.1008303 0.02028824 -0.03678615 -0.05092357 0.02339115 0.0240892 0.01927424 -0.04905638 0.04213179 -0.0005581815 -0.1033197 0.1318134 -0.09368698 0.08528423 -0.02169268 -0.03422522 -0.04499843 0.08258557 -0.03179607 -0.01507485 0.02660707 -0.02762987 0.0711119 -0.08979131 -0.008203965## $test
## [1] "Ljung-Box test"
##
## $K
## [1] 48
##
## $chi.square
## [1] 143.4359
##
## $df
## [1] 48
##
## $pval
## [1] 0.00000000001924017acf(e$res)
pacf(e$res)
dev.off()## null device
## 1#Compare Spectral Densities
sims = 5
SpecDen = parzen.wge(us$positivity, plot = "FALSE")
plot(SpecDen$freq, SpecDen$pzgram, type = "l", lwd = 6)
for( i in 1: sims)
{
SpecDen2 = parzen.wge(gen.aruma.wge(246,s = 7, phi = e$phi, plot ="FALSE"), plot = "FALSE")
lines(SpecDen2$freq,SpecDen2$pzgram, lwd = 2, col = "red")
}
#Compare ACFs
sims = 5
ACF = acf(us$positivity, plot = "FALSE")
plot(ACF$lag ,ACF$acf , type = "l", lwd = 6)
for( i in 1: sims)
{
ACF2 = acf(gen.aruma.wge(246, s = 7, phi = e$phi, plot = "FALSE"), plot = "FALSE")
lines(ACF2$lag ,ACF2$acf, lwd = 2, col = "red")
}
#Compare Generated Realizations
eGen = gen.aruma.wge(246, s = 7, phi = e$phi, vara = e$avar)
plotts.sample.wge(eGen)## $autplt
## [1] 1.0000000 0.9848470 0.9714933 0.9750294 0.9660213 0.9445656 0.9409767
## [8] 0.9406086 0.9164590 0.8947200 0.8887335 0.8694487 0.8382897 0.8244560
## [15] 0.8136619 0.7805097 0.7512321 0.7387271 0.7140631 0.6782452 0.6590084
## [22] 0.6421241 0.6032096 0.5685146 0.5501356 0.5211394
##
## $freq
## [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
## [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
## [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
## [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
## [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
## [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
## [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
## [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
## [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
## [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
## [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
## [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
## [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
## [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
## [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
## [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
## [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
##
## $db
## [1] 17.4883412 17.8759163 3.3192846 -3.8525491 -0.3197144 -21.6714749
## [7] -22.8063713 -10.2392453 -12.0673705 -12.2214139 -13.9805703 -18.7730531
## [13] -18.8355978 -14.2310658 -18.3550959 -24.3257081 -18.1713112 -18.9508838
## [19] -21.0905082 -22.5490074 -22.9798514 -27.1738870 -22.4637158 -26.0337988
## [25] -16.3206084 -8.5546634 -15.3803215 -19.2883534 -17.8332777 -18.1899535
## [31] -26.4043277 -23.8962853 -22.7948280 -16.8336539 -6.1759350 -17.3886878
## [37] -26.3712238 -30.6911752 -30.0931877 -25.9624947 -23.5436142 -25.9359153
## [43] -25.3700146 -23.9993226 -25.8558264 -26.4139571 -23.7900585 -30.4384324
## [49] -30.0297904 -27.3820589 -30.1075893 -19.9753587 -21.7443303 -19.2880151
## [55] -26.0065533 -34.3012107 -21.6232692 -24.7615997 -27.7416393 -26.4778644
## [61] -26.6516604 -21.3567198 -22.4913621 -21.6969354 -20.1932738 -20.4355328
## [67] -17.8545927 -16.6305613 -12.4553796 -2.0051496 -6.0649279 -18.2027733
## [73] -22.1637145 -21.1998955 -24.6167114 -27.6646930 -30.7627963 -28.9387210
## [79] -37.0841490 -35.0228031 -33.7229017 -30.6743821 -35.4618060 -32.0714755
## [85] -30.0881307 -29.4898690 -37.9232247 -40.4322393 -36.4403348 -37.8417324
## [91] -31.3156666 -36.6701417 -32.4208118 -37.4804681 -31.6757217 -28.1692244
## [97] -30.6926889 -37.2504701 -43.4538675 -29.2992604 -34.3024720 -31.7441604
## [103] -28.5361267 -22.1853685 -24.1843615 -18.4839851 -27.5289666 -33.4879999
## [109] -39.8646199 -33.4828829 -34.6946852 -38.9555566 -35.0056370 -46.1802585
## [115] -38.4655865 -34.6454560 -39.9518903 -36.9388807 -24.9050466 -41.8672693
## [121] -39.9673409 -38.8068328 -34.7290668
##
## $dbz
## [1] 13.1505985 12.8508351 12.3476953 11.6357495 10.7071467 9.5513936
## [7] 8.1551522 6.5022292 4.5742086 2.3529001 -0.1723734 -2.9851629
## [13] -5.9943454 -8.9461610 -11.3785406 -12.8764690 -13.4514878 -13.3828185
## [19] -12.9305630 -12.3235211 -11.7616150 -11.3810556 -11.2483531 -11.3764193
## [25] -11.7386223 -12.2740811 -12.8899611 -13.4736343 -13.9256735 -14.2014856
## [31] -14.3245535 -14.3576413 -14.3642960 -14.3912229 -14.4710836 -14.6318058
## [37] -14.9027654 -15.3154018 -15.8998317 -16.6793579 -17.6628684 -18.8327983
## [43] -20.1258665 -21.4114139 -22.4978380 -23.2156270 -23.5363512 -23.5726581
## [49] -23.4591612 -23.2742952 -23.0485995 -22.8036367 -22.5760745 -22.4147948
## [55] -22.3601135 -22.4148600 -22.5064242 -22.4471161 -21.9635668 -20.8921253
## [61] -19.3571507 -17.6410485 -15.9700130 -14.4625165 -13.1665084 -12.0964210
## [67] -11.2533126 -10.6339936 -10.2348609 -10.0534452 -10.0890256 -10.3428890
## [73] -10.8184658 -11.5214051 -12.4595315 -13.6424492 -15.0801763 -16.7792701
## [79] -18.7327507 -20.8959192 -23.1374554 -25.1828816 -26.6789790 -27.4915621
## [85] -27.8444184 -28.0356068 -28.2163791 -28.4034520 -28.5572936 -28.6450395
## [91] -28.6696151 -28.6616936 -28.6512341 -28.6425116 -28.6051411 -28.4842864
## [97] -28.2288276 -27.8241888 -27.3050207 -26.7379019 -26.1926286 -25.7235865
## [103] -25.3653040 -25.1357367 -25.0412690 -25.0807136 -25.2475968 -25.5307047
## [109] -25.9130085 -26.3693131 -26.8636964 -27.3492391 -27.7738304 -28.0941293
## [115] -28.2925476 -28.3849520 -28.4114379 -28.4166945 -28.4327401 -28.4708364
## [121] -28.5224704 -28.5672362 -28.5849455plotts.sample.wge(us$positivity)## $autplt
## [1] 1.0000000 0.9470164 0.9341278 0.9260111 0.9198625 0.9088272 0.8909389
## [8] 0.8836047 0.8738064 0.8521265 0.8325090 0.8001230 0.7724089 0.7494712
## [15] 0.7330396 0.6973509 0.6643364 0.6399202 0.6205127 0.5900024 0.5527729
## [22] 0.5218086 0.4844797 0.4575403 0.4134570 0.3899807
##
## $freq
## [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
## [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
## [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
## [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
## [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
## [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
## [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
## [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
## [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
## [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
## [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
## [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
## [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
## [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
## [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
## [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
## [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
##
## $db
## [1] 16.928178 17.041622 10.954209 3.357603 1.069112 -3.238369
## [7] -11.427364 -14.831916 -13.181967 -26.879328 -10.640049 -13.104348
## [13] -6.897898 -10.769553 -12.587656 -16.736482 -15.880304 -19.355323
## [19] -19.864320 -18.877313 -12.131779 -10.974583 -16.180375 -7.015352
## [25] -8.700809 -8.412062 -7.545175 -6.229323 -11.252124 -16.095149
## [31] -15.372565 -27.552336 -10.346212 -9.057794 -11.461366 -7.109201
## [37] -14.240281 -13.075078 -23.877820 -15.330952 -14.475847 -20.765059
## [43] -13.129351 -25.189770 -19.078147 -17.562843 -19.295449 -16.258703
## [49] -12.990975 -10.740763 -14.551855 -9.308655 -7.444421 -5.100658
## [55] -8.399431 -9.881040 -17.975943 -17.462588 -17.722631 -13.792923
## [61] -13.293199 -9.316456 -15.869391 -14.012814 -14.465016 -12.635251
## [67] -8.756403 -15.569710 -11.708858 -22.389341 -7.883153 -18.081602
## [73] -16.077320 -13.828997 -15.999068 -10.429352 -26.201640 -19.135996
## [79] -14.542536 -16.213531 -16.985225 -20.728879 -15.016682 -31.829738
## [85] -14.128733 -13.936455 -8.979219 -13.464535 -15.769438 -12.144105
## [91] -14.497446 -17.272748 -12.667286 -11.379894 -13.465174 -9.903788
## [97] -13.261966 -12.491846 -11.986603 -11.921460 -13.694377 -12.698553
## [103] -16.097000 -15.677901 -11.088504 -11.527014 -9.506795 -12.681243
## [109] -14.748674 -18.413266 -16.493124 -18.659561 -18.956181 -17.539240
## [115] -23.088025 -19.681462 -12.465460 -11.043980 -10.130884 -8.698403
## [121] -17.478083 -17.827845 -23.241743
##
## $dbz
## [1] 12.8975742 12.6104343 12.1288216 11.4480751 10.5614946 9.4601976
## [7] 8.1330155 6.5665777 4.7459620 2.6568532 0.2916097 -2.3347667
## [13] -5.1453437 -7.9215132 -10.2286515 -11.5956090 -11.9705143 -11.6983430
## [19] -11.1158574 -10.4346224 -9.7877324 -9.2573489 -8.8845067 -8.6785829
## [25] -8.6271903 -8.7045894 -8.8785406 -9.1161191 -9.3884024 -9.6733156
## [31] -9.9562335 -10.2290872 -10.4895332 -10.7412048 -10.9945308 -11.2664674
## [37] -11.5777577 -11.9475163 -12.3857167 -12.8839971 -13.4053244 -13.8761561
## [43] -14.1925522 -14.2554748 -14.0264019 -13.5536343 -12.9396293 -12.2881747
## [49] -11.6760017 -11.1508180 -10.7394607 -10.4558243 -10.3057268 -10.2889439
## [55] -10.3991649 -10.6224555 -10.9347537 -11.2993911 -11.6668145 -11.9799663
## [61] -12.1877268 -12.2626791 -12.2120173 -12.0727873 -11.8955891 -11.7286297
## [67] -11.6092436 -11.5622931 -11.6019261 -11.7339670 -11.9576927 -12.2666738
## [73] -12.6487555 -13.0854284 -13.5510287 -14.0125543 -14.4313384 -14.7679001
## [79] -14.9901220 -15.0821781 -15.0493866 -14.9156322 -14.7148182 -14.4811756
## [85] -14.2423645 -14.0164221 -13.8117305 -13.6288608 -13.4634734 -13.3096700
## [91] -13.1631127 -13.0231056 -12.8930232 -12.7790349 -12.6877269 -12.6235629
## [97] -12.5870636 -12.5742955 -12.5779435 -12.5898586 -12.6044544 -12.6217762
## [103] -12.6489228 -12.6990700 -12.7883489 -12.9315227 -13.1373601 -13.4041478
## [109] -13.7156159 -14.0383220 -14.3232947 -14.5159238 -14.5746988 -14.4899695
## [115] -14.2884211 -14.0193071 -13.7340281 -13.4719928 -13.2561370 -13.0948841
## [121] -12.9868227 -12.9259627 -12.9064607# Check performance
preds <- fore.aruma.wge(us$positivity, phi = e$phi, theta = e$theta, s=7, n.ahead = 7, lastn = T, limits = F)
eval_model(us$positivity,preds$f, model_name = 'ARMA(18,1) With Weekly/Monthly Trend', AIC_val = e$aic) #ASE .000178
preds <- fore.aruma.wge(us$positivity, phi = e$phi, theta = e$theta, s=7, n.ahead = 90, lastn = T, limits = F)
eval_model(us$positivity,preds$f, model_name = 'ARMA(18,1) With Weekly/Monthly Trend', AIC_val = e$aic) #ASE .000707
rolling <- rolling_ASE(us, e, s=7, horizon=7, model_name = 'ARMA(18,1) with Weekly Trend') #ASE .000093## [1] "test window: 240:246, train window: 210:239"## [1] "Window ASE: 0.000160550780936603"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 210:216, train window: 180:209"## [1] "Window ASE: 0.0000856798636025404"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 180:186, train window: 150:179"## [1] "Window ASE: 0.0000230595258989943"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 150:156, train window: 120:149"## [1] "Window ASE: 0.000230459645617526"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 120:126, train window: 90:119"## [1] "Window ASE: 0.000032800419583193"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 90:96, train window: 60:89"## [1] "Window ASE: 0.0000252534592547995"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## Warning: Removed 89 row(s) containing missing values (geom_path).Signal Plus Noise
The signal plus noise model performs very well in the 7 day window, but very poorly over the 90 day window. This is a deterministic signal with a stationary mean model, and our results aren’t surprising since it doesn’t account for the weekly trend well.
preds <- fore.sigplusnoise.wge(us$positivity, max.p = 12, n.ahead = 7, limits=F)##
## Coefficients of Original polynomial:
## 0.4763 0.1833 0.1282 0.1498 0.0814 -0.0857 0.0905 0.1286 -0.0125 0.0792 -0.1415 -0.1573
##
## Factor Roots Abs Recip System Freq
## 1-1.9128B+0.9191B^2 1.0406+-0.0722i 0.9587 0.0110
## 1-1.3370B+0.8117B^2 0.8236+-0.7441i 0.9009 0.1169
## 1-0.3817B+0.7741B^2 0.2465+-1.1096i 0.8798 0.2152
## 1+1.3262B+0.7343B^2 -0.9031+-0.7392i 0.8569 0.3908
## 1+0.4065B+0.7162B^2 -0.2838+-1.1470i 0.8463 0.2886
## 1+1.4225B+0.5177B^2 -1.3737+-0.2106i 0.7195 0.4758
##
## 
eval_model(us$positivity,preds$f, preds$ul, preds$ll,'SigPlusNoise',AIC_val = 0) #ASE = .000169
preds <- fore.sigplusnoise.wge(us$positivity, max.p = 12, n.ahead = 90, limits=F)##
## Coefficients of Original polynomial:
## 0.4899 0.1825 0.1263 0.1449 0.0783 -0.0882 0.0842 0.1163 -0.0191 0.0749 -0.1362 -0.1288
##
## Factor Roots Abs Recip System Freq
## 1-1.9009B+0.9072B^2 1.0476+-0.0687i 0.9525 0.0104
## 1-1.3188B+0.7894B^2 0.8353+-0.7543i 0.8885 0.1169
## 1-0.3724B+0.7559B^2 0.2463+-1.1235i 0.8694 0.2157
## 1+1.3275B+0.7233B^2 -0.9177+-0.7351i 0.8505 0.3925
## 1+0.4049B+0.6962B^2 -0.2908+-1.1627i 0.8344 0.2890
## 1+1.3697B+0.4724B^2 -1.4497+-0.1229i 0.6873 0.4865
##
## 
eval_model(us$positivity,preds$f, preds$ul, preds$ll,'SigPlusNoise',AIC_val = 0) #ASE = .002678
rolling <- rolling_ASE(us, e, s=7, horizon=7, model_name = 'SigPlusNoise', model_type = 'SigPlusNoise', p = 15) #ASE .00012## [1] "test window: 240:246, train window: 210:239"
##
## Coefficients of Original polynomial:
## 0.2342 -0.4855 -0.1223 -0.2171 -0.0969 -0.1119 0.3218
##
## Factor Roots Abs Recip System Freq
## 1-1.1778B+0.9139B^2 0.6444+-0.8240i 0.9560 0.1444
## 1+0.2335B+0.8023B^2 -0.1455+-1.1069i 0.8957 0.2708
## 1+1.3852B+0.6500B^2 -1.0656+-0.6349i 0.8062 0.4145
## 1-0.6751B 1.4812 0.6751 0.0000
##
## 
## [1] "Window ASE: 0.000208585888080764"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 210:216, train window: 180:209"
## [1] "Window ASE: 0.0000342786896513529"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 180:186, train window: 150:179"
##
## Coefficients of Original polynomial:
## -0.1222 -0.3150
##
## Factor Roots Abs Recip System Freq
## 1+0.1222B+0.3150B^2 -0.1939+-1.7712i 0.5612 0.2674
##
## 
## [1] "Window ASE: 0.0000312625829956592"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 150:156, train window: 120:149"
## [1] "Window ASE: 0.000268642564641623"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 120:126, train window: 90:119"
## [1] "Window ASE: 0.000057908435916914"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 90:96, train window: 60:89"
##
## Coefficients of Original polynomial:
## -0.1703 -0.2635 -0.1903 0.1825 0.0224 0.1150 0.1611 0.0453 0.0858 -0.2458 -0.4970 -0.0954 -0.3683
##
## Factor Roots Abs Recip System Freq
## 1-1.8801B+0.9434B^2 0.9965+-0.2590i 0.9713 0.0405
## 1-1.2359B+0.9403B^2 0.6572+-0.7947i 0.9697 0.1400
## 1+0.9641B -1.0373 0.9641 0.5000
## 1+1.6980B+0.9042B^2 -0.9389+-0.4737i 0.9509 0.4256
## 1+0.9115B+0.8715B^2 -0.5229+-0.9349i 0.9336 0.3312
## 1-0.0174B+0.7523B^2 0.0116+-1.1529i 0.8674 0.2484
## 1-0.2699B+0.7264B^2 0.1857+-1.1585i 0.8523 0.2247
##
## 
## [1] "Window ASE: 0.000120756392201658"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## Warning: Removed 89 row(s) containing missing values (geom_path).
Neural Network
The neural network model performs fairly well in the short term, and has the best rolling window ASE out of the bunch using minimal parameter tuning (setting lags = 7, and doing 5 fold cross validation). It is, however, the most difficult to explain, using Multilayer Perceptrons (MLP) for prediction. The example below uses two hidden nodes, both of which are autoregression lags, which apply a weigh parameter to the training input, and to each other, resulting in the prediction outcome.
ts_us <- ts(us$positivity[1:239], start = '1')
x = mlp(ts_us, lags = 7, hd.auto.type = 'cv')## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).plot(x)
preds <- predict(x, 7)
plot(preds)
eval_model(us$positivity,preds$mean, model_name = 'NNFOR', AIC_val = 0) #ASE = .00021## Warning: attributes are not identical across measure variables;
## they will be dropped
ts_us <- ts(us$positivity[1:156], start = '1')
x = mlp(ts_us, hd.auto.type = 'cv', sel.lag = T)
plot(x)
preds <- predict(x, 90)
plot(preds)
eval_model(us$positivity,preds$mean,model_name = 'NNFOR', AIC_val = 0) #ASE = .002771## Warning: attributes are not identical across measure variables;
## they will be dropped
rolling <- rolling_ASE(us, e, s=0, horizon=7, model_name = 'NNFOR', model_type = 'NNFOR') #ASE .00009## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.0000362940126580657"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 210:216, train window: 180:209"## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).
## [1] "Window ASE: 0.0000153520775693192"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 180:186, train window: 150:179"## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).
## [1] "Window ASE: 0.0000429492839377769"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 150:156, train window: 120:149"## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).
## [1] "Window ASE: 0.000239811808026264"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 120:126, train window: 90:119"## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).
## [1] "Window ASE: 0.0000669272791843012"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 90:96, train window: 60:89"## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).
## [1] "Window ASE: 0.000169772044069679"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## Warning: Removed 89 row(s) containing missing values (geom_path).
Multivariate Analysis
Taking what we learned from our univariate analysis, let’s see if we can add some explanatory variables to explain more difference in the data. For Louisiana, I have selected average temperature, wind speed, precipitation, as well as significant events that I think would result in increased gatherings. These include holidays, holiday weekends, the two weeks of school start dates across Louisiana, protests resulting from the Black Lives Matter movement, and various voting days.
For the US dataset, I will only be using the significant dates data, as the weather data only makes sense when applied to a specific region, such as Louisiana.
Louisiana
MLP
la_XDF <- data.frame(
temp = ts(xdf$AvgTempF),
humidity = ts(xdf$AvgHumidityPercent),
wind = ts(xdf$AvgWindSpeedMPH),
precip <- ts(xdf$PrecipitationIN),
gathering <- ts(xdf$GatheringBinary)
)
## 7 day forecast
ts_la <- ts(la$positivity[1:239], start = '1')
x = mlp(ts_la, hd.auto.type = 'cv', lags=7, xreg = la_XDF)
plot(x)
preds <- forecast(x, h=7, xreg=la_XDF)
plot(preds)
eval_model(la$positivity, preds$mean, model_name = 'NNFOR', AIC_val = 0) #ASE = .000664## Warning: attributes are not identical across measure variables;
## they will be dropped
## 90 day forecast
ts_la <- ts(la$positivity[1:156], start = '1')
x = mlp(ts_la, hd.auto.type = 'cv', lags=7, xreg = la_XDF)
plot(x)
preds <- forecast(x, h=90, xreg=la_XDF)
plot(preds)
eval_model(la$positivity, preds$mean, model_name = 'NNFOR', AIC_val = 0) #ASE = .000176## Warning: attributes are not identical across measure variables;
## they will be dropped
rolling <- rolling_ASE(la, s=0, horizon=7, model_name = 'NNFOR', model_type = 'NNFOR') #ASE .000535## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.000971578708755374"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 210:216, train window: 180:209"
## [1] "Window ASE: 0.000125575980051967"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 180:186, train window: 150:179"
## [1] "Window ASE: 0.000156638653876302"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 150:156, train window: 120:149"
## [1] "Window ASE: 0.00015584632792789"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 120:126, train window: 90:119"
## [1] "Window ASE: 0.00150668585334667"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 90:96, train window: 60:89"
## [1] "Window ASE: 0.000545042330189"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## Warning: Removed 89 row(s) containing missing values (geom_path).
VAR
VARselect(la$positivity[1:239], lag.max = 7, type = "both", season = 7, exogen = la_XDF[1:239,]) #AIC = -8.5337900068## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 2 2 1 2
##
## $criteria
## 1 2 3 4 5
## AIC(n) -8.5264373480 -8.5337900068 -8.5269742401 -8.5200371542 -8.5124988698
## HQ(n) -8.4425560122 -8.4439171471 -8.4311098564 -8.4181812465 -8.4046514381
## SC(n) -8.3184445755 -8.3109406078 -8.2892682145 -8.2674745019 -8.2450795910
## FPE(n) 0.0001981888 0.0001967436 0.0001980968 0.0001994845 0.0002010039
## 6 7
## AIC(n) -8.5091143397 -8.5027303904
## HQ(n) -8.3952753840 -8.3828999107
## SC(n) -8.2268384342 -8.2055978583
## FPE(n) 0.0002016965 0.0002030007#VARselect picks p=2 using AIC and p=1 using BIC
vfit=VAR(cbind(Positivity = la$positivity, la_XDF)[1:239,], p=1,type='both', season = 7)
preds=predict(vfit,n.ahead=7)
eval_model(la$positivity, preds$fcst$Positivity[,1], pred_ul = preds$fcst$Positivity[,3],
pred_ll = preds$fcst$Positivity[,2], model_name = 'VAR', AIC_val = -8.5337900068) #ASE = .000650
vfit=VAR(cbind(Positivity = la$positivity, la_XDF)[1:156,], p=1,type='both', season = 7)
preds=predict(vfit,n.ahead=90)
eval_model(la$positivity, preds$fcst$Positivity[,1], pred_ul = preds$fcst$Positivity[,3],
pred_ll = preds$fcst$Positivity[,2], model_name = 'VAR', AIC_val = -8.5337900068) #ASE = .000625
rolling <- rolling_ASE(la, s=7, horizon=7, model_name = 'VAR', model_type = 'VAR', p = 1, df_XDF = la_XDF) #ASE .000281## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.000469447094936274"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 210:216, train window: 180:209"
## [1] "Window ASE: 0.000150722025793936"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 180:186, train window: 150:179"
## [1] "Window ASE: 0.0000558290163363635"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 150:156, train window: 120:149"
## [1] "Window ASE: 0.000134574853250694"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 120:126, train window: 90:119"
## [1] "Window ASE: 0.000486713927394731"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 90:96, train window: 60:89"
## [1] "Window ASE: 0.000389727851279653"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## Warning: Removed 89 row(s) containing missing values (geom_path).
United States
MLP
us_XDF <- data.frame(
gathering <- ts(xdf$GatheringBinary)
)
ts_us <- ts(us$positivity[1:239], start = '1')
x = mlp(ts_us, lags = 7, hd.auto.type = 'cv', xreg = us_XDF)## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).plot(x)
preds <- predict(x, 7)
plot(preds)
eval_model(us$positivity,preds$mean,model_name='NNFOR', AIC_val = 0) #ASE = .000204## Warning: attributes are not identical across measure variables;
## they will be dropped
ts_us <- ts(us$positivity[1:156], start = '1')
x = mlp(ts_us, hd.auto.type = 'cv', lags = 7, xreg = us_XDF)## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).plot(x)
preds <- predict(x, 90)
plot(preds)
eval_model(us$positivity,preds$mean, model_name='NNFOR', AIC_val = 0) #ASE = .004178## Warning: attributes are not identical across measure variables;
## they will be dropped
rolling <- rolling_ASE(us, s=0, horizon=7, model_name = 'NNFOR', model_type = 'NNFOR') #ASE .000098## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.0000365472078758298"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 210:216, train window: 180:209"## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).
## [1] "Window ASE: 0.0000153293647402558"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 180:186, train window: 150:179"## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).
## [1] "Window ASE: 0.0000439697478267936"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 150:156, train window: 120:149"## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).
## [1] "Window ASE: 0.000256916030066162"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 120:126, train window: 90:119"## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).
## [1] "Window ASE: 0.0000639204210989819"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length## [1] "test window: 90:96, train window: 60:89"## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).
## [1] "Window ASE: 0.000140189117737396"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## Warning: Removed 89 row(s) containing missing values (geom_path).
VAR
VARselect(us$positivity[1:239], lag.max = 7, type = "both", season = 7, exogen = us_XDF[1:239,]) #AIC = -8.5191652146## Warning in VARselect(us$positivity[1:239], lag.max = 7, type = "both", season = 7, : No column names supplied in exogen, using: exo1 , instead.## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 6 4 4 6
##
## $criteria
## 1 2 3 4 5
## AIC(n) -8.3318103491 -8.4588431215 -8.4875954850 -8.5030175995 -8.4956045603
## HQ(n) -8.2718951092 -8.3929363577 -8.4156971972 -8.4251277877 -8.4117232245
## SC(n) -8.1832440830 -8.2954202288 -8.3093159658 -8.3098814537 -8.2876117879
## FPE(n) 0.0002407487 0.0002120323 0.0002060271 0.0002028791 0.0002043947
## 6 7
## AIC(n) -8.5035667501 -8.4950952290
## HQ(n) -8.4136938903 -8.3992308453
## SC(n) -8.2807173511 -8.2573892034
## FPE(n) 0.0002027806 0.0002045136#VARselect picks p=4 using BIC and p=6 using AIC... p=6 returns NA values, so using p=4
vfit=VAR(cbind(Positivity = us$positivity, la_XDF)[1:239,], p=4,type='both', season = 7)
preds=predict(vfit,n.ahead=7)
eval_model(la$positivity, preds$fcst$Positivity[,1], pred_ul = preds$fcst$Positivity[,3],
pred_ll = preds$fcst$Positivity[,2], model_name = 'VAR', AIC_val = -8.5337900068) #ASE = .001218
vfit=VAR(cbind(Positivity = us$positivity, la_XDF)[1:156,], p=4,type='both', season = 7)
preds=predict(vfit,n.ahead=90)
eval_model(us$positivity, preds$fcst$Positivity[,1], pred_ul = preds$fcst$Positivity[,3],
pred_ll = preds$fcst$Positivity[,2], model_name = 'VAR', AIC_val = -8.5337900068) #ASE = .001590
rolling <- rolling_ASE(us, s=7, horizon=7, model_name = 'VAR', model_type = 'VAR', p = 4, df_XDF = us_XDF) #ASE .000122## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.000170145130090431"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 210:216, train window: 180:209"
## [1] "Window ASE: 0.0000490511908879048"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 180:186, train window: 150:179"
## [1] "Window ASE: 0.0000488062219779536"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 150:156, train window: 120:149"
## [1] "Window ASE: 0.000299162276131377"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 120:126, train window: 90:119"
## [1] "Window ASE: 0.0000907668571318613"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 90:96, train window: 60:89"
## [1] "Window ASE: 0.0000714752472492872"## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## Warning: Removed 89 row(s) containing missing values (geom_path).
Model Comparison
| Type | Model | Rolling ASE | Dataset |
|---|---|---|---|
| Univariate | Louisiana | AR(15) Weekly | .000277 |
| Univariate | Louisiana | ARIMA(14,0,1) Weekly | .000679 |
| Univariate | Louisiana | SigPlusNoise | .000418 |
| Univariate | Louisiana | MLP | .000637 |
| Multivariate | Louisiana | MLP | .000601 |
| Multivariate | Louisiana | VAR | .000281 |
| — | — | — | — |
| Univariate | National | AR(9,1) | .000096 |
| Univariate | National | ARIMA(18,0,1) Weekly | .000093 |
| Univariate | National | SigPlusNoise | .000120 |
| Univariate | National | MLP | .000090 |
| Multivariate | National | MLP | .000098 |
| Multivariate | National | VAR | .000122 |
Final Model Selection
LA model: Weekly AR(15)
𝜑(1−.682𝐵−.23〖3𝐵〗2−.179𝐵3+.144𝐵4+.035𝐵5−.073𝐵6+.644𝐵7−.364𝐵8−.082𝐵9+.066𝐵10+.128𝐵11−.015𝐵12−.216𝐵13+.341𝐵14−.306𝐵15)(1−𝐵7)=𝑎𝑡 - 7-day ASE = . 000647 - 90-day ASE = . 000202 - Rolling ASE = .000277 - AIC = -8.488985
US model: AR(9,1)
𝜑(1−1.271𝐵+.163𝐵2+.070𝐵3−.066𝐵4−.003𝐵5+.164𝐵6−.161𝐵7−.114𝐵8+.229𝐵9)=(1−.882𝐵)𝑎𝑡 - 7-day ASE = .000121 - 90-day ASE = .001150 - Rolling ASE = .000096 - AIC = –8.625600
Summary
In both cases, the Univariate models outperformed their multivariate counterparts. This may suggest that I selected poorly correlated features to include in the model. But as a testament to the time series algorithms, they perform just as well as the much more complicated neural networks.
e <- est.arma.wge(la.s7, p=15, q=0)##
## Coefficients of Original polynomial:
## 0.6819 0.2329 0.1795 -0.1438 0.0353 0.0725 -0.6438 0.3639 0.0825 -0.0664 -0.1285 0.0148 0.2163 -0.3415 0.3063
##
## Factor Roots Abs Recip System Freq
## 1-1.8995B+0.9868B^2 0.9625+-0.2950i 0.9934 0.0473
## 1+0.9661B+0.9242B^2 -0.5226+-0.8994i 0.9614 0.3338
## 1+1.3916B+0.9115B^2 -0.7634+-0.7172i 0.9547 0.3800
## 1+1.8666B+0.8970B^2 -1.0405+-0.1795i 0.9471 0.4728
## 1-0.9455B 1.0576 0.9455 0.0000
## 1-0.1935B+0.8186B^2 0.1182+-1.0989i 0.9048 0.2329
## 1-1.3104B+0.7391B^2 0.8865+-0.7531i 0.8597 0.1121
## 1-0.5571B+0.7181B^2 0.3879+-1.1145i 0.8474 0.1967
##
## preds <- fore.aruma.wge(la$positivity, phi = e$phi, theta = e$theta, s=7, n.ahead = 30, lastn = F, limits = F)
ggplot() +
geom_line(aes(us$date, us$positivity), color='#F8766D') +
geom_line(aes(
seq(as.Date("2020/11/12"), by = "day", length.out = 30),
preds$f), color='#00BFC4') + labs(title = 'Lousiana 30 day predictions', y = 'Positivity Rate', x='Date')
e <- est.arma.wge(us$positivity, p=9, q=1)##
## Coefficients of Original polynomial:
## 1.2711 -0.1631 -0.0699 0.0662 0.0028 -0.1636 0.1612 0.1144 -0.2295
##
## Factor Roots Abs Recip System Freq
## 1-1.9861B+0.9886B^2 1.0045+-0.0501i 0.9943 0.0079
## 1-0.0326B+0.7214B^2 0.0226+-1.1772i 0.8493 0.2469
## 1-1.1925B+0.6943B^2 0.8588+-0.8384i 0.8332 0.1231
## 1+1.2068B+0.6319B^2 -0.9548+-0.8190i 0.7949 0.3872
## 1+0.7333B -1.3636 0.7333 0.5000
##
## preds <- fore.aruma.wge(us$positivity, phi = e$phi, theta = e$theta, s=7, n.ahead = 30, lastn = F, limits = F)
ggplot() +
geom_line(aes(us$date, us$positivity), color='#F8766D') +
geom_line(aes(
seq(as.Date("2020/11/12"), by = "day", length.out = 30),
preds$f), color='#00BFC4')+ labs(title = 'US 30 day predictions', y = 'Positivity Rate', x='Date')